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1.  INTRODUCTION 


1-1 .  General  Aspects 

Naval  Electronics  systems  process  extremely  large  quantities  of  data  today  and 
it  is  expected  that  demands  on  high  speed  real-time  signal  processing  will 
increase  substantially  in  the  next  decade.  Requisite  device  speed  and  densi¬ 
ties  demand  that  critical  feature  size  be  reduced  to  submicron  or 
ultrasubmicron  scales. 

This  reduction  places  uncommon  demands  on  the  quality  of  the  materials  used  in 
device  fabrication.  For  example,  in  present  bipolar  devices,  crystallographic 
microdefects,  precipitates,  clusters,  and  deep  level  traps  can  result  in  micro¬ 
plasma  and  junction  leakage,  current  channeling,  nonuniforra  heating,  degraded 
pn  junction  breakdown,  etc.  Taken  in  combination  these  introduce  limits  on 
size  and/or  chip  yield.  In  silicon  MOS  devices,  similar  defects  lead  to  gener¬ 
ation-recombination  centers  and  leakage  current  sources.  Here  refresh  times 
in  dynamic  RAMs,  threshold  shifts  in  static  RAMs ,  maximum  operating  tempera¬ 
tures,  etc.,  are  degraded.  The  situation  in  gallium  arsenide  is  even  more 
difficult.  Defects  in  the  substrate  material  are  severely  limiting  yield, 
raising  costs,  and  causing  time  scales  for  projected  commercial  and  military 
applications  to  be  pushed  back  several  years.  This  is  true  for  both  GaAs 
digital  devices  and  the  new  MMIC  devices. 

Material  variations  arise  during  all  stages  of  device  fabrication,  beginning 
with  crystal  growth.  In  most  methods  of  synthetic  production  of  single  crys¬ 
tals,  the  crystal  grows  slowly  from  a  fluid  nutrient.  Being  a  fluid,  various 
hydrodynamic  mechanisms  drive  it's  motion.  The  fluid  motion  is  of  concern  to 
crystal  growers  because  they  determine  the  transport  of  impurities,  heat,  etc., 
to  the  grown  crystal,  and  the  stresses  on  the  crystal. 

The  present  study  involves  the  application  of  numerical  techniques  to  the  most 
popular  method  of  growing  crystals:  the  Czochralski  (Cz)  process.  The  Cz 
method  is  the  principle  method  by  which  silicon  and  gallium  arsenide  boules 
are  obtained  in  the  microelectronics  industry.  The  study  was  supported  by  the 
Office  of  Naval  Research  through  the  DESAT  program.  The  motivation  for  the 
program  is  clear.  It  is  currently  accepted  that  the  properties  of  the  melt 
are  the  primary  factors  affecting  the  quality  of  silicon  and  gallium  arsenide 
substrates.  The  properties  of  the  melt  are  themselves  determined  by  a  broad 
and  interacting  matrix  of  variables  that  include  crystal  rotation,  crucible 
rotation,  pull  rate,  crucible  temperature,  pressure,  orientation,  etc.  These 
input  variables  have  been  traditionally  operator  controlled,  and  the  quality 
of  the  boule  was  often  dependent  upon  the  past  experience  of  the  operator,  and 
resulting  trial  and  error  procedures.  Because  of  the  high  demands  of  both  the 
commercial  and  military  markets,  as  well  as  international  competition,  the 
control  of  the  growth  variables  and  the  manner  in  which  a  variation  of  one 
affects  another,  systematic  techniques  are  needed  for  optimization  of  crystal 
growth.  The  very  first  step  in  this  optimization  is  the  development  of  numer¬ 
ical  codes  that  closely  replicate  the  growth  procedures.  The  second  is  the 
development  of  software  that  is  coupled  to  the  actual  growth  apparatus  and  is 
capable  of  being  self -correcting.  At  this  point  in  time,  both  of  these  proce¬ 
dures  are  under  develop  ent.  But  the  degree  to  which  this  coupled  approach  is 
useful  is  dependent  upon  the  constraints  of  the  equations  governing  the  fluid 
flow. 
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Pragmatically,  present  crystal  growers  require  a  three  dimensional  transient 
algorithm  that  describes  crystal  growth  and  is  capable  of  handling  such 
external  effects  such  as  the  presence  of  magnetic  field,  electric  fields, 
encapsulation,  melt  replenishment,  etc.  Such  algorithms  have  been  developed 
at  Scientific  Research  Associates,  Inc.  (SRA)  and  have  been  applied  to  silicon 
and  gallium  arsenide . 

The  equations  solved  in  the  SRA  study  are  the  continuity,  momentum  and  energy 
balance  equations  for  the  melt,  as  well  as  the  species  concentration 
equation.  For  the  crystal,  the  energy  balance  equation  is  solved.  The 
governing  equations  are  formulated  using  primitive  variables  and  are  not 
thereby  restricted  to  limitations  of  two  dimensional  symmetry.  The  crystal 
and  melt  phases  are  coupled  through  thermal  and  dynamic  conditions  applied 
along  the  melt-crystal  and  melt-free  surface  boundaries.  For  most  of  the 
calculations  the  shape  of  the  interface  is  determined  through  solutions  to 
differential  equations  that  govern  the  kinematic  conditions  at  the  specific 
interface,  and  are  not  assumed  a  priori.  The  numerical  procedures  are 
obtained  by  a  well  documented  procedure  known  as  the  consistently  split 
Linearized  Block  Implicit  (LBI)  scheme,  originally  developed  at  SRA.  While 
SRA's  numerical  studies  of  crystal  growth  are  perhaps  the  most  general,  to 
date,  the  extent  to  which  the  procedures  and  results  are  of  relevance  should 
be  placed  in  perspective  at  this  early  stage  in  the  discussion. 

All  numerical  studies  of  Cz  growth  involve  simplifying  assumptions  and  result 
in  placing  the  studies  in  different  groups,  or  categories  as  summarized  in 
Table  1.1.1.  The  first  groups  involves  only  studies  of  the  thermal -hydrody¬ 
namics  of  the  melt  phase.  That  is,  the  energy  and  momentum  balance  equations 
for  the  melt  were  solved.  The  assumptions  for  this  group  are  generally:  1) 
planar  free  surface;  2)  planar  melt-crystal  interface,  and;  3)  uniform 
temperature  distribution  in  the  crystal  phase.  The  second  group  involves  the 
studies  of  temperature  distribution  in  the  melt  phase  and/or  the  crystal 
phase.  The  influence  of  hydrodynamics  of  the  melt  on  the  temperature  is 
assumed  to  be  small;  and  only  the  energy  balance  equation  is  solved.  However, 
attention  in  this  group  was  paid  to  the  shape  of  the  free  surface  as  well  as 
that  of  the  melt  surface  interface.  The  third  group  includes  that  of  the 
first  but  also  includes  energy  balance  within  the  crystal.  The  melt  and 
crystal  phases  are  coupled  together  via  thermal  conditions  at  the  melt -crystal 
interface.  The  analysis  on  this  group  is  difficult,  and  the  shape  of  the  free 
surface  and  the  melt  surface  interface  is  not  considered. 

The  present  work  extends  that  of  the  last  group,  plus:  i)  the  shape  of  the  free 
surface  and  includes  the  mensicus  at  the  crystal-melt-argon  contact  point,  ii) 
the  shape  of  the  melt-crystal  interface,  iii)  normalized  dopant  distribution 
in  the  melt,  iv)  dopant  distribution  near  the  melt-crystal  interface,  v)  the 
effect  of  replenishment,  vi)  the  dependence  of  melt  properties  on  crucible 
size,  vii)  the  dependence  of  melt  properties  on  crystal  size,  viii)  the 
dependence  of  melt  properties  on  crystal  penetration  of  the  melt. 

The  structure  of  the  paper  is  as  follows:  Section  1.2  presents  a  brief 
discussion  of  CZ  growth.  This  is  followed  by  a  historical  summary  of  the 
earlier  studies.  Section  2  describes  the  mathematical  formulation  of  the  CZ 
growth  problem.  In  Section  3  numerical  details  are  discussed,  including  a 
description  of  the  body-fitted  coordinate  algorithm,  transient  capability  and 
CRAY  vertorization.  Section  4  contains  a  description  of  model  'numerically 


based'  calculations.  Application  to  silicon  is  contained  in  Sections  5  and  6, 
with  the  latter  devoted  to  three  dimensional  studies.  Conclusions  and 
Recommendations  of  future  directions  is  contained  in  Section  7. 

1-2.  Fundamentals  of  Czochralski  growth 

The  basic  principle  of  growing  a  single  crystal  by  the  Czochralski  method  is 
shown  in  Figure  1.1.1.  The  crucible  is  initially  filled  with  the  material 
e.g.,  polycrystalline  Si,  to  be  crystallized.  Thermal  energy  is  supplied  by  a 
heater  surrounding  the  crucible  thereby  raising  the  temperature  of  the  solid 
above  its  melting  temperature.  A  seed  crystal,  at  the  end  of  a  rod,  is  then 
dipped  into  the  melt  and,  after  appropriate  start-up  procedures,  slowly 
withdraw  from  the  melt. 

Under  suitable  temperature  conditions  recrystallizaton  in  the  form  of  single 
crystals  occurs.  This  method,  bears  Czochralski 's  name,  who  in  1916  [1], 
established  it  as  a  means  to  determine  the  crystallization  velocity  of 
metals.  For  the  growth  of  silicon  (Si)  and  germanium  (Ge)  single  crystals, 
the  typical  pulling  rate  is  less  than  3.5mm/min.  For  a  2  inch  diameter 
gallium  arsenide  crystal  is  about  0.15  mm/min. 

To  keep  the  feeding  material  molten,  the  crucible  is  maintained  at  a 
temperature  above  the  melting  point.  The  resulting  axial  and  radial 
temperature  gradients  within  the  melt  cause  natural  convection.  In  order  to 
improve  uniformity,  the  growing  crystal  is  rotated  as  it  is  pulled.  The 
rotating  crystal  acts  as  a  centrifugal  pump  sucking-up  the  melt  fluid  axially, 
while  ejecting  it  tangentially. 

As  the  crystal  rotates,  a  viscous  shear  layer  (Ekman  layer)  is  formed  under  the 
crystal,  tending  to  isolate  the  growing  solid/melt  interface  from  the  melt  near 
the  bottom  of  the  crucible.  The  radially  outward  flow,  due  to  the  rotating 
crystal,  is  met  by  a  countering  radial  inflow  fluid  driven  by  the  hot  crucible 
wall . 

A  resultant  downflow  occurs  at  some  radial  distance,  whose  position  is 
determined  by  the  relative  strength  of  the  two  rotational  flows  and  the 
buoyancy- driven  convection. 

Rotation  of  the  crucible,  which  is  also  sometimes  utilized,  is  implemented  to 
smooth  out  additional  thermal  asymmetries  that  arise  from  irregularities  in 
heating.  This  also  introduces  a  centrifugal  flow.  The  combined  effects  of 
crystal  and  crucible  rotation,  even  in  an  iso-thermal  fluid,  are  complex  and 
depend  on  whether  the  two  are  iso-rotating  or  counter-rotating. 

With  iso-rotation,  a  Taylor -Proudman  cell  can  be  established  beneath  the 
crystal,  with  a  detached  shear  layer  separating  it  from  the  outer  region.  The 
outer  region  effectively  rotates  as  a  solid  body  with  the  crucible  [2]. 
Counter-rotation  yields  a  more  complex  flow  pattern  beneath  the  crystal  with  a 
stagnation  layer  separating  the  flows  driven  by  the  crystal  from  the  flows 
driven  by  crucible  rotation  [2,3],  Typical  crystal  and  crucible  rotation 
rates  are  in  the  range  of  1  to  50  and  0.5  to  30rpm,  respectively. 


In  addition  to  the  natural  and  forced  convection,  there  is  thermo-capillary 
flow  near  the  free  surface.  Since  the  surface  tension  coefficients  of  most 


Cz  melts  (for  Si,  it  is  -0.43  dyne/cm/k)  vary  with  temperature,  a  non-zero 
temperature  gradient  along  the  free  surface  will  cause  a  nonuniform  surface 
tension  force,  resulting  in  a  finite  shear  stress  at  the  free  surface.  For 
example,  for  a  melt  with  negative  coefficient  of  surface  tension,  the  surface 
tension  diminishes  with  increasing  radius,  provided  the  crucible  temperature 
is  higher  than  the  melt  temperature.  The  resultant  effect  is  that  the  melt 
beneath  the  free  surface  is  pulled  toward  the  crystal.  Motion  induced  in  this 
manner  is  called  Marangoni  flow,  and  is  further  complicated  by  the  nonplanar 
nature  of  the  free  surface,  the  loss  of  heat  from  the  free  surface  and  the 
consequent  destabilizing  vertical  temperature  gradient  that  gives  rise  to  a 
Benard-like  convection. 


In  addition,  the  shape  of  the  melt-crystal  interface  is  controlled  by  both  the 
coupled  effects  of  heat  transfer  throughout  the  crystal  growing  system  and  the 
shape  of  the  free  surface.  Capillarity  and  the  wetting  of  the  crystal  by  its 
melt  determine  the  position  and  shape  of  the  meniscus. 

1- 3 .  Background 

Recent  tutorial  and  review  articles  [4-8]  display  the  nature  of  such 
convective  motions  and  their  influence  on  the  quality  of  the  crystal . 

Intensive  research  to  understand  the  hydrodynamics  of  the  Cz  growth  system 
have  been  carried  out  in  the  past  decade.  As  discussed  earlier,  numerical 
investigation  of  Cz  growth  can  be  classified  into  three  groups  (see  Table 
1.1.1).  Group  A  involves  the  studies  on  the  hydrodynamics  of  the  melt  phase 
with  the  assumption  of  a  planar  free  surface  and  melt  crystal  interface. 

Group  B  involves  studies  of  temperature  distribution  within  both  melt  and 
crystal  phases.  Group  C  studies  combine  both  that  of  Groups  A  and  B,  without 
the  assumption  of  planar  interfaces. 

The  Group  A  assumptions  include:  1)  The  melt-crystal  interface  and  free 

surface  are  planar;  2)  The  crucible  and  crystal  are  at  uniform  temperatures 

and  rotating  symmetrically;  3)  Crystal  phase  is  not  studied;  4)The  pulling 
velocity  is  zero.  The  first  assumption  ignores  the  kinetics  of  the  surfaces 
which  are  important  when  the  solute  distribution  and  thermal  gradients  near 
the  interface  are  to  be  determined.  The  second  assumption  rests  with 
formulation  limitations  as  it  is  restricts  to  axisymmetric  situations.  Since 
the  crystal  phase  is  not  studied  in  this  group,  effects  due  to  heat  loss  at 

the  vertical  surface  of  the  crystal  cannot  be  incorporated  at  the  planar 

melt-crystal  interface. 

The  pioneering  work  on  numerical  simulation  of  convective  flow  was  carried  out 
by  Kobayashi  [9,11]  and  Langlois  [10]  using  a  model  that  employed  a  vorticity- 
stream  function  procedure.  In  these  models,  the  axial  and  radial  velocities 
were  replaced  by  a  stream  function  through  the  use  of  the  continuity 
equation.  The  pressure  was  explicitly  eliminated  by  introducing  the  vorticity 
function,  and  the  governing  equations  were  then  written  in  finite  difference 
form.  Kobayashi  [9,11]  further  assumed  a  steady  state  and  solved  each 
equation  successively.  Since  the  unknowns  were  coupled,  the  entire  procedure 
was  repeated  iteratively  until  prespecified  convergence  criteria  was  met. 
Langlois  [10]  formulated  the  problem  as  a  transient  and  used  a  forward  time 
explicit  scheme.  The  stream  function  (continuity)  was  calculated  at  each  time 
step  by  a  successive  over-relaxation  method.  This  numerical  method  (either 
steady  or  unsteady)  was  then  used  to  simulate  convective  flows  in  Cz  systems 


(also  in  a  floating  zone  system,  e.g.,  [12])  with  various  nondimens ional 
parameters  [13-16].  Instead  of  using  a  finite  difference  and  SOR  methods  for 
the  vorticity-stream  function  formulation,  Crochet  and  Wouters  [17]  used 
finite  element  methods  to  study  convective  flows  due  to  various  rotational 
rates  and  buoyancy  forces.  The  numerical  method  has  been  demonstrated  to  be 
able  to  handle  oscillatory  natural  convective  flow.  These  papers  have  made 
significant  contributions  to  the  understanding  of  the  hydrodynamics  of  the 
melt . 

In  the  Group  B  studies,  only  the  temperature  distribution  in  the  melt  and/or 
crystal  phase  is  solved.  A  study  by  Williams  and  Reusser  [18]  used  a  finite 
element  method  to  determine  the  temperature  distribution  in  the  crystal. 

Their  work  assumed  a  planar  melt-crystal  interface  and  specified  the  boundary 
conditions  from  the  measured  surface  temperatures.  This  approach  is  not 
general  and  is  unsuited  to  predictive  purposes.  More  recently,  Derby  et  al., 

[19]  included  a  meniscus  and  nonplanar  interfaces  in  two  dimensional,  steady 
state,  small  scale  prototype  calculations  applicable  to  liquid  encapsulated 
Czochralski  (LEC)  growth.  The  thermal  equation  in  each  phase  was  solved  using 
a  finite  element  technique.  However,  phenomena  associated  with  fluid  flow 
were  totally  ignored.  Thus  the  shape  of  the  meniscus  was  simplified  to  the 
well  known  Young-Laplace  equation,  for  which  an  analytical  solution  exists 

[20]  .  The  isotherm  at  the  melt  temperature  was  used  to  define  the 
melt-crystal  interface.  Since  the  flow  field  was  not  solved,  a  solute 
distribution  cannot  be  determined  from  this  analysis.  A  study  by  Ramachandran 
and  Dudukovic  [21]  used  an  iterative  finite  element  method  to  predict  the 
temperature  distribution  and  the  position  of  the  melt-crystal  interface  in  the 
crystal  phase.  The  method  accounted  for  the  radiative  heat  transfer  between 
the  various  'surfaces'  presented  to  the  crystal  pulling  apparatus. 

In  Group  C  studies,  Chang  and  Brown  [21]  formulated  the  vertical  Bridgeman 
system  using  primitive  variables  and  approximated  the  solution  using  a 
Galerkin  finite  element  method.  In  their  study,  the  steady  state  two 
dimensional  Navier-Stokes  equations  are  solved  for  the  case  in  which  flow  is 
due  to  free  convection  only.  The  isotherm  at  the  melt  temperature  is 
considered  as  the  melt-crystal  interface  and  is  approximated  by  a  set  of 
Hermite  cubic  polynominals .  The  coefficients  associated  with  the  field 
variables  and  interface  shape  formed  a  nonlinear  set  of  algebraic  equations. 

In  the  vertical  Bridgeman  system,  the  melt  fluid  is  continuously  fed  from  the 
bottom  while  the  crystal  is  formed  at  the  top.  No  free  surface  exists  in  this 
configuration.  The  range  of  parameters  studied  in  these  groups  is  summarized 
in  Table  1.1.2. 

Because  of  the  difficulty  in  making  measurements  in  the  melt  (e.g.  opaqueness, 
high  temperature  of  the  melt),  simulation  experiments  [23-27]  have  be">n 
performed  to  study  melt  behavior  under  Cz  conditions.  The  crystal  is 
simulated  by  a  rotating  disc,  and  the  crucible  by  a  heated  cylindrical 
container.  The  working  fluid  (simulating  the  melt)  is  taken  to  be  water 
[23-23,  27]  and  NaN03  [26]  so  that  the  Prandtl  number  is  of  the  same  order 
as  the  oxide  melts.  Photographs  from  Miller  [25]  show  a  'wheel  spoke'  pattern 
at  the  free  surface.  The  'ripples'  are  believed  to  be  due  to  Marangoni 
convection  resulting  from  surface  tension  gradients  along  the  free  surface. 

At  high  crystal  rotation  rates,  three  distinct  zones  are  observed  [25]:  the 
stagnant  layer  at  the  bottom,  the  thermally  driven  zone  in  the  middle,  and  the 
crystal  driven  zone  near  the  crystal.  Lamprecht  et  al . ,  [26]  put  tracer 
particles  in  NaN03 ,  so  that  when  the  fluid  was  illuminated  with  a  sheet  of 
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light,  photographs  could  be  taken.  By  measuring  the  relative  positions  of  the 
tracer's  particles  between  time  frames,  radial  velocities  were  measured.  The 
temperature  distribution  was  measured  with  thermocouples.  Marangoni 
convection  was  not  observed  because  the  surface  tension  coefficient  for 
NaN03  is  small  (-0.07  dyne/cm/k  at  300°C,  -0.17  dyne/cm/k  for  water  at 
50°C).  In  another  experiment  done  with  water  (Jone  [27]),  an  array  of 
thermocouples  was  located  1  cm  below  a  rotating  disc  (simulated  crystal)  to 
measure  the  mean  and  fluctuating  temperature  when  the  rate  of  rotation  of  the 
disc  was  increased.  Observed  experimental  flow  patterns  and  temperature 
variations  qualitatively  agree  with  the  melt  behavior  under  Czochralski 
conditions . 

Another  aspect  of  crystal  growth  being  studied  is  how  the  dopant  impurities 
and  heat  are  transported  to  the  melt-crystal  interface  as  a  result  of  the 
convective  flows.  The  flows  due  to  the  rotating  crystal  and  the  resultant 
solute  distribution  produced  by  segregation  at  the  melt  crystal  interface  have 
been  treated  by  Wilson  and  others  [28-32]  using  the  similarity  solutions  of 
Von  Karman  [33]  for  the  fluid  flow.  The  melt-crystal  interface  is  prescribed 
to  be  a  planar  disc  normal  to  its  axis  of  rotation.  A  suitable  coordinate  can 
be  employed  to  accommodate  the  nonplanar  crystal  surface.  The  importance  of 
these  studies  is  that  it  is  the  behavior  of  the  fluid  adjacent  to  the  growing 
crystal  that  is  the  most  influential  in  controlling  the  solute  uptake.  The 
bulk  of  the  melt  provides  the  reservoir  of  solute.  The  spatial  distribution 
of  the  solute  in  the  melt  is  assumed  to  be  uniform.  Originated  from  Burton  et 
al.,  [34]  an  analytical  formula  can  be  written  for  the  effective  segregation 
coefficient  as  a  function  of  equilibrium  segregation  coefficient  and  a 
parameter,  A.  The  latter  can  be  interpreted  as  a  concentration  boundary 
layer  thickness.  Several  numerical  correlations  for  the  parameter,  A,  as  a 
function  of  the  rotational  rate  and  pulling  rate  of  the  crystal  were  developed 
in  [28,32].  Additionally,  Wilson  and  Favier  [30,31,33]  also  included  a  time 
dependent  sinuodial  varying  pulling  rate  in  their  model  to  study  compositional 
inhomogenities  in  the  crystal.  The  rotating  disk  model  [28-32,34]  was  limited 
to  cases  with  a  uniform  spatial  impurity  distribution.  Molecular  oxygen  which 
comes  from  a  continuous  corrosion  of  the  crucible  walls,  is  transported 
throughout  the  melt  and  the  majority  of  it  evaporates  from  the  free  surface 
with  some  fraction  incorporated  in  the  grown  crystal.  In  this  case,  the 
spatial  distribution  within  the  melt  is  not  uniform  and  thus  this  model  is  not 
applicable. 

As  discussed  earlier  the  goal  of  the  present  study  was  to  overcome  many  of  the 
limitations  associated  with  earlier  studies  and  to  develop  numerical 
procedures  that  are  capable  of  accurately  predicting  for  Cz  processes:  1)  the 
thermal  hydralic  behavior  of  the  melt,  2)  the  thermal  behavior  of  the  crystal; 
3)  the  shapes  of  the  free  surface  and  melt-crystal  interface;  and  4)  the 
species  distribution  in  the  melt  fluid. 


2.  MATHEMATICAL  ANALYSIS 


2-1.  Introduction 


In  tht  present  study,  the  governing  equations  for  the  liquid  state  are  the 
conservation  equations  for  mass,  momentum  and  energy  of  the  melt.  The  only 
equation  considered  for  the  crystal  is  the  heat  conduction  equation.  The  melt 
and  crystal  phases  are  coupled  through  appropriate  thermal  conditions  applied 
along  the  interfaces.  The  interfacial  positions  are  governed  by  nonlinear 
kinematic  conditions,  as  discussed  below.  The  solutions  to  the  melt  flow 
field  and  the  temperature  distribution  inside  the  grown  crystal,  along  with 
the  positions  of  the  interface,  are  obtained  by  a  well  documented  numerical 
procedure  known  as  the  consistently  split  linearized  block  implicit  (LBI) 
scheme  [35]  (see  Appendix  A). 

The  addition  of  free  surface  boundary  conditions  introduces  nonlinear  effects, 
since  the  location  of  free  surface  is  not  known  a  priori,  and  must  emerge  as 
part  of  the  solution.  The  unknown  surface  location  is  obtained  from  a  partial 
differential  equation  boundary  condition  that  is  coupled  to  flow  within  the 
melt.  At  the  melt-solid  interface,  the  conservation  equations  for  the  melt 
are  coupled  to  the  thermal  conservation  equation  for  the  crystal  through  the 
energy  balance  condition. 

In  the  calculations  of  the  species  concentration  in  the  melt,  species  are 
assumed  to  be  transported  by  two  mechanisms;  i)  melt  fluid  motion,  and  ii) 
species  diffusion.  Concentration  gradients  induced  by  the  temperature  rise  at 
the  melt-crystal  interface  due  to  both  latent  heat  of  crystallization  and  the 
curvature  of  the  melt-crystal  interface  are  assumed  to  be  small.  Thus  motion 
due  to  these  concentration  gradients  is  assumed  to  be  small  and  is  neglected. 

2.2.  Governing  Equations 

The  equations  describing  flow  within  the  melt  and  energy  transport  within  the 
crystal  are  discussed  below.  Each  of  these  equations  are  expressed,  in  two 
forms;  a)  in  dimensional  form,  and  b)  in  dimensionless  form. 

The  first  equation  of  interest  is  that  governing  the  distribution  of 
temperature  within  the  solid  crystal.  The  equation  for  heat  transport  is: 
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In  the  above  Ts  is  the  temperature  of  the  crystal,  up'  the  speed  at  which 
the  crystal  is  pulled,  as  is  the  thermal  diffusivity  which  is  defined  as 
ks/^gcps>  and  ks  is  the  thermal  conductivity,  cps  is  the  specific  heat 
capacity  and  ps  is  the  density.  The  dimensionless  form  of  the  equations 
is  clear.  0S  represents  T/Tm  where  Tm  is  the  temperature  of  the 
melt.  up  -  up ' /Ug .  The  reference  velocity  is  us  -  Rs  u>s .  Rs  is  the  radius 


of  the  crystal.  ws  is  the  angular  frequency  of  rotation  of  the  crystal. 

The  dimensionless  number  Pe  is  referred  to  as  the  Peclet  number  and  is  given 
by  Pe  -  usrs/as.  Note  that  the  crystal  radius  is  taken  as  the  reference 
length,  and  the  reference  time  is  l/«s. 

The  governing  equations  for  the  melt  are  as  follows:  (1)  First  there  is  the 
continuity  equation 


p  +  V  .  (  pljf  )  =  0 
-2-p  +  V  .  (  pu  )  =0 


where  p'  is  the  density  of  the  melt  and  u'  its  velocity.  In  dimensionless 
form  the  density  is  normalized  to  the  density  of  the  melt  at  the  melt 
temperature.  The  second  melt  equation  is  the  momentum  balance  equation,  or 
more  commonly  the  Navier- Stokes  equation 
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Here  F'  is  a  force  per  unit  mass  and  P'  is  the  pressure,  u  is  the  kinematic 
viscosity  which  is  defined  as  p/p' .  In  dimensionless  units,  the  reference 
pressure  is  pus2.  The  dimensionless  Reynolds  number  Re  is  given  by 
pusRs,  where  p  is  the  viscosity. 

The  third  melt  equation  is  the  energy  equation 
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In  the  above  a£  is  the  thermal  diffusivity.  In  dimensionless  form,  Pr  is 
the  Prandtl  number,  and  is  given  by  »VQj.  The  fourth  melt  equation, 
species  transport  equation,  is 
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where  Sc  is  referred  to  as  the  Schmidt  number,  and  is  given  by  p/D,  and  D  is 
species  diffusivity. 


The  force  per  unit  mass  term  in  equation  (3)  is  given  by 

F'  =  -fl  (  1  +  M  T  -  T  )  ) 
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where  g  is  the  acceleration  due  to  gravity  and  /9  is  the  volumetric  expansion 
coefficient  of  the  melt.  Fr  is  the  Froude  number  and  defined  as  us2/gRs, 
while  Gr  is  the  Grashof  number  and  is  defined  as  g/3ATR®/i/2 .  Here,  the 

Bousinesq  approximation  is  adopted  to  account  for  the  effects  of  natural 
convection.  The  formulation  of  the  crystal  growth  problem  allows  for  the 
incorporation  of  compressible  fluid  properties.  While  not  likely  to  be 
important  for  a  melt  that  is  exposed  to  low  pressure,  significant  departures 
may  occur  under  either  high  pressure  and/or  encapsulation.  The  incorporation 
of  compressibility  into  the  algorithm  requires  the  presence  of  an  equation  of 
state.  Care  must  be  taken  in  the  choice  of  such  an  equation,  in  that  it  must 
yield  results  that  are  virtually  identical  to  those  obtained  by  others  in  the 
incompressible  limit.  The  equation  of  state  used  below  has  been  successfully 
applied  to  compressible  fluids,  and  yields  the  same  results  as  others  in  the 
incompressible  limit  [36],  This  equation  is: 


(7) 


where  R  is  the  universal  gas  constant. 

Further,  according  to  the  analysis  given  in  [36],  if  Tr  were  the  local  total 
temperature  in  the  fluid,  equation  (7)  would  represent  the  perfect  gas  law  for 
an  ideal  fluid.  However,  if  Tr  is  taken  as  a  fictictious  constant  temperature 
of  high  enough  value  so  that  the  Mach  number  Ma  based  upon  u^  and  Tr  is 
small ,  i . e . , 


Ma 


l*|l 


/ 


«  1.0 


Y  R(Tr  -  |  / 2C 


(8) 


'<1 


ft 


then,  equations  (2)  and  (3)  will  represent  the  governing  equations  for  nearly 
incompressible  flows.  Since,  in  the  present  work,  the  melt  is  considered  as 
nearly  incompressible,  all  the  cases  investigated  in  the  present  work  have  Mach 
numbers  on  the  order  of  10'2. 


Equations  (1)  through  (5)  are  written  in  vector  form.  The  specific  dependent 
variables  to  be  determined  are  the  temperatures,  pressure,  concentration  and 
velocity  components  in  the  direction  of  the  cylindrical  polar  coordinate 
lines.  In  addition,  the  locations  of  the  free  surface  and  melt-crystal 
interface  will  be  determined  along  with  the  solution  of  the  above  equations .  An 
important  and  novel  feature  of  the  study  is  that  a  time -dependent,  generalized 
coordinate  transformation  is  introduced  to  accommodate  the  motion  and  the  shape 
of  the  free  surface  and  melt-crystal  interface. 


2-3.  Free  Surface  Boundary  Conditions 


The  free  surface  defined  here  is  the  interface  between  the  melt  and  a  gas.  In 
this  situation,  the  stress  tangential  to  the  surface  is  approximated  as  zero. 
The  stress  normal  to  the  surface  must  exactly  balance  any  externally  applied 
normal  stress.  With  surface  tension  included,  the  stress  conditions  are 
derived  from, 
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where  nj  are  the  components  of  the  unit  vector  normal  to  the  surface,  tjj  is 
the  stress  tensor  and  y  is  the  force  due  to  surface  tension.  For  a  Newtonian 
fluid  and  axisymmetric  flow,  the  stress  conditions  can  be  simplified  as  [38]: 


(i)  Tangential  stress  condition  (along  the  azimuthal  direction) 
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(ii)  Tangential  stress  condition  (along  the  radial  direction) 
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(iii)  Normal  stress  condition 
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where  Bo  is  the  Bond  number  and  is  defined  as  pgRs/7,  and  determines  the 
ratio  of  the  surface  tension  to  gravational  forces,  there  r)  is  the  free 
surface  elevation  and  r)'  ,  rj' '  are  first  and  second  derivatives  with  respect 
to  r,  respectively.  nt ,  n2 ,  and  n3  are  r,  p ,  z  components  of  a  unit 
normal  vector  at  the  free  surface.  The  stress  conditions  are  also  called 
dynamic  conditions  because  they  are  the  force  balance  conditions  which  must  be 
satisfied  along  the  free  surface  at  all  times. 

In  addition  to  the  dynamic  conditions,  a  kinematic  constraint  is  used  to 
determine  the  position  of  the  free  surface.  This  condition  is  based  upon  the 
constraint  that  the  normal  velocity  of  the  fluid  at  the  interface  is  the  same 
as  the  normal  velocity  of  the  free  surface.  In  Eulerian  form,  it  can  be 
expressed  as, 

^■(n-z)  =  o  (13) 

where  D_  is  the  substantial  deriviative  [44] .  The  location  of  the  free 
DT 

interface  is  determined  by  time  integration  of  the  above  equation. 


2-4 .  Velocity  Boundary  Conditions 


Velocity  boundary  conditions  along  the  melt-crystal  interface  and  crucible 
wall  are  also  required.  Along  the  melt-crystal  interface,  'no-slip'  and  mass 
balance  conditions  are  imposed.  The  'no-slip'  condition  implies 
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and  the  mass  balance  condition  is, 
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where  r ^ ,  n^  are  the  unit  tangential  and  normal  vectors  at  the 
melt-solid  interface.  The  velocity  boundary  conditions  along  the  crucible 
wall  are 
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where  r c ,  nc  are  the  unit  tangential  and  normal  vectors  at  crucible 
wall,  e.g.,  the  tangential  vector  at  the  crucible  bottom  wall  is  a  horizontal 
vector,  while  the  tangential  vector  at  the  crucible  side  wall  is  a  vertical 
vector . 


A  special  kind  of  velocity  boundary  condition  needs  to  be  used  at  the  junction 
where  the  moving  surface  meets  the  stationary  wall.  Physically,  'no-slip' 
conditions  should  be  used  at  these  points;  however,  it  will  cause  numerical 


problems  if  u  is  zero  at  the  junction.  Here,  according  to  the  kinematic 
condition,  the  surface  will  never  move  at  these  points.  Experience  teaches 
that  the  surface  should  move  according  to  the  rotating  boundary  conditions. 
Thus,  at  points  where  the  moving  points  meet  the  crucible  side  wall,  we  allow 
velocities  tangential  to  the  walls  to  slip,  and  velocities  normal  to  the  wall 
to  be  'no  slip' .  The  slip  velocities  are  obtained  from  the  kinematic 
condition.  At  every  time  step,  the  kinematic  conditions  are  applied  at  the 
surface,  except  at  the  junction.  The  location  of  the  surface  at  the  wall  is 
obtained  from  a  linear  extrapolation  from  the  interior  points. 

2-5.  Melt-crvstal  Interface  Boundary  Conditions 

There  are  two  conditions  that  need  to  be  satisfied  along  the  melt-crystal 
interface;  1)  the  temperature  must  be  equal  to  the  equilibrium  melting 
temperature,  and  2)  the  net  heat  flux  across  the  interface  must  equal  to  the 
latent  heat  absorbed  per  unit  volume  during  solidification.  The  first 
condition  states  that  the  melt-solid  interface  is  a  constant  temperature 
surface  i.e.,  T  -  T  .  The  second  condition  can  be  written  into  the  form  of  a 
kinematic  condition  and  used  to  determine  the  position  of  the  melt-crystal 
interface.  We  will  now  develop  the  second  condition. 

The  heat,  AQ| ,  delivered  to  or  removed  away  from  the  melt -crystal 
interface  by  phase  i  in  a  time  interval,  At,  can  be  written  as 

AQ  .“An  At 
i  n  n.i 


where  qn  ^  is  the  hear  flux  normal  tc  the  interface  and  An  is  the  cross- 
sectional  area  perpendicular  to  the  heat  flux.  Using  Fourier’s  law  for  the 
heat  flux,  the  net  heat  delivered  to  the  melt-crystal  interface  is 
Q  -  AQ^ - AQS ,  i.e.. 
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where  n  is  the  unit  normal  vector  at  the  melt  crystal  interface.  Here,  heat 
transported  through  radiation  is  ignored  This  net  heat,  according  to  the  law 
of  conservation  of  energy,  should  equal  to  the  latent  heat  gained  due  to 
crystallization  process  which  is 
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where  A r)*  is  the  change  of  interfacial  positions  in  time,  At;  and  hs^  is 
the  latent  of  crystallization 


Coupling  the  net  heat  flux  across  the  interface 
yields : 


to  the  latent  heat  absorbed 
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To  treat  equation  (19)  specifically,  we  note  that  foi  the  Cz  problem,  the 


crystal  is  pulled  out  Up  at  a  constant  velocity  (within  a  reasonable  time 
period),  i.e.,  the  melt-crystal  interface  is  rising  at  a  constant  rate.  If 

rj*  is  written  as  A rj  -  Afj  where  A»j  is  the  area  averaged  change  of  interfa¬ 
cial  positions  and  Arj  accounts  for  the  spatial  variations  of  change  of  the 
interfacial  positions,  then  Aq/At'  -  Arj*/At'  +  A»j/A t' .  Further,  if  At'  is 
small,  the  above  quantities  can  be  approximated  by  time  derivatives,  then  the 
above  equation  can  be  written  as 
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where  dji  —  up'  is  the  crystal  pulling  velocity.  In  nondimens ional  form, 
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where  Pe  and  S  are  Peclet  and  Stefan  numbers,  respectively.  S  is  defined  as 
hsj>/CpSAT  and  k  is  the  ratio  of  melt  thermal  conductivity.  Note  that 
in  crystal  growth  problems,  Up  is  the  crystal  pulling  rate  which  is  exter¬ 
nally  controlled;  thus,  dri/d t  is  a  measure  of  nonuniform  crystal  growth  (in 
both  space  and  time).  Furthermore,  the  shape  of  the  melt-crystal  interface  is 
of  great  interest  in  crystal  growth  process  control  because  it  affects  the 
distribution  of  dopants  and  unwanted  impurities  in  the  crystal,  and  is 
responsible  for  residual  strains  in  the  crystal. 


In  the  present  model,  the  crystal  pulling  velocity  is  determined  from 
integrating  equation  (21),  i.e., 
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Physically,  the  crystal  pulling  velocity  is  obtained  from  the  vc lume  averaged 
net  heat  across  the  interface.  Since  we  have  assumed  a  constant  crvstal 
diameter  growth  in  a  steady  state  situation,  the  total  net  heat  removed  from 
the  melt-crystal  structure  will  be  from  the  heat  gained  due  to  crystalliza¬ 
tion.  In  steady  state  calculations,  this  relation  is  exact,  but  in  transient 
calculations,  we  have  assumed  the  volume  average  of  the  non-uniform  growth  is 
small . 

2 - 6 .  Thermal  Boundary  Conditions 

Three  types  of  thermal  boundary  conditions  are  involved  in  our  analysis  for 
the  liquid  and  solid  phases.  They  are,  1)  constant  temperature  boundary 


1J 


condition,  2)  constant  heat  flux  boundary  condition,  and  3)  constant  heat 
radiation  boundary  condition.  For  the  liquid  phase,  constant  temperature 
boundary  condition  are  always  applied  at  the  crucible  walls  and  crystal 
surface.  A  combination  of  heat  flux  and  heat  radiation  boundary  conditions 
can  be  applied  at  the  free  surface.  For  the  solid  phase,  constant  temperature 
boundary  condition  is  applied  at  the  melt  crystal  interface.  These  boundary 
conditions  are  outlined  next. 

1)  Constant  temperature  boundary  condition 

Liquid:  9g  -  at  crystal  surface 

-  9C  at  crucible  wall  (23) 

Solid:  9S  -  0m  at  melt  solid  interface 


2)  Heat  flux  condition 


Since  argon  gas  is  blowing  along  the  crystal  surface  and  free  surface, 
convective  heat  transfer  may  be  important,  thus: 


Liquid: 

ii 

.  kg,Nu(6,-ea) 

at  free  surface 

Solid: 

39  _ 

at  crystal  surface 

_ _ s  _ 

3n 
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(above  free  surface) 

(24) 


kgi  is  the  ratio  of  conductivity  of  the  gas  to  conductivity  of  Phase  i 
(j>  or  s) ,  Nu  is  the  Nusselt  number  which  is  defined  as  hRs/kg  where  h 
is  the  heat  transfer  coefficient  of  the  gas  and  k  is  the  thermal  conductivity 
of  the  gas.  0a  is  the  ambient  temperature.  The  Nusselt  number  can  be 
obtained  through  some  empirical  relations,  e.g.  in  Appendix  B.  Note,  at  the 
top  of  crystal,  a  zero  flux  condition  is  used  to  simulate  an  infinitely  long 
grown  crystal. 

3)  Heat  Radiation  Condition 

In  a  crystal  growing  system  involving  N  surfaces  with  different  emissivities 
at  different  temperatures,  the  mutual  radiation  interchange  is  adequately 
treated  by  introducing  the  concept  of  geometrical  shape  factor,  f^,- .  In  an 
ideal  situation,  f^j  is  1.  The  net  heat  flux  radiated  from  the  kth  surface 
is  expressed  by, 
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When  considering  only  the  direct  radiation  interchange,  a  is  the  Stefan- 
Boltzmann  constant,  Tj  and  ej  are  the  temperature  and  emissivity  of  the 
jth  surface.  For  the  present  application,  we  only  consider  heat  radiation 
from  the  free  surface  to  the  gas  (argon)  or  heat  radiation  from  the  crystal 
surface  above  the  free  surface  to  the  gas  phase,  and  assume  a  constant  ambient 
temperature  above  the  free  surface.  A  better  approximation  for  the  ambient 


temperature  can  be  used  if  the  heat  transfer  and  flow  conditions  in  argon  gas 
are  known.  The  heat  flux  radiated  can  be  simplified  as, 
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Nondimens ionalizing  the  above  equation  and  introducing  a  Biot  number,  H,  the 
above  equation  becomes 
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where  Rc  is  the  nondimens ional  crucible  radius,  and  H  is  defined  as 

eoRsTm3Ai. 


2-7 .  Species  Boundary  Conditions 

At  the  melt-crystal  interface,  the  flux  condition  on  the  concentration  is, 
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where  is  the  equilibrium  distribution  coefficient  for  species  k.  The 

suction  parameter,  a,  is  defined  as  Up//(wv) .  We  have  assumed  that  diffu¬ 
sion  into  the  crystal  is  negligible.  For  most  species,  in  Si  growth,  the 
coefficient  aReV2  Sc  is  on  the  order  of  100  with  Sc  on  the  order  of  50. 
Numerically,  it  means  a  very  fine  grid  must  be  used  near  the  interface  to 
resolve  the  boundary  layer.  Physically,  the  boundary  layer  for  the  convec • 
tive-diffussive  species  equation  is  approximately  Sc^/^  times  thinner  than 
the  momentum  boundary  layer.  To  reduce  computational  efforts,  we  do  not 
resolve  the  boundary  layer  of  the  species  equation  in  the  bulk  flow 
calculations.  Instead,  we  use  the  zero  flux  condition: 
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To  obtain  the  resolution  of  the  species  distribution  near  the  interface,  a 
solute  transport  submodel  has  been  developed.  This  submodel  makes  use  of  the 
solution  obtained  from  the  'bulk  flow'  calculations. 


At  the  crucible  walls,  as  well  as  the  free  surface,  a  zero  flux  condition  is 
also  used.  Here  we  assume  the  dissolution  and  vaporization  rates  at  the 
crucible  walls  and  free  surface  respectively,  are  small.  In  order  to 
calculate  the  oxygen  concentration  in  the  melt,  a  model  for  the  rate  of 
dissolution  and  vaporization  of  oxygen  needs  to  be  developed  (either  from 
experimental  measurements  or  derived  from  crucible  properties  and  argon  gas 
flow  dynamics).  This  is  beyond  the  scope  of  the  present  effort. 


2  -  8 .  Species  Transport  Submodel 

It  was  found  [28-32,34]  that  the  behavior  of  the  fluid  adjacent  to  the  growing 
crystal  interface  is  a  primary  factor  in  controlling  the  species  uptake  to  the 
crystal.  In  our  model,  we  assume  that  species  diffusion  and  forced  convection 
(due  to  rotation  of  the  crystal)  are  important  in  a  thin  region  with  thickness 
A  near  the  melt-crystal  interface.  Solute  impurities  in  the  melt  and  from 
tfie  crucible  wall  are  transported  by  "bulk"  flow  and  reach  the  boundary  of 
this  thin  region.  The  governing  equations  (2)  and  (3) ,  can  be  simplified  to  a 
set  of  ordinary  differential  equations  [33]  to  provide  a  flow  field  for  the 
species  equation  (5)  in  this  region. 

In  the  thin  region,  if  the  natural  convection  is  ignored,  the  nondimens ional- 
ized  axisymmetric  continuity  [2]  and  Navier-Stokes  equations  [3]  can  be 

(29) 
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(31) 


(32) 


where  equation  (29)  is  the  continuity  equation.  Here  equation  (30)  througn 
(32)  are  the  three  components  of  Navier-Stokes  equation  and  u  -  rwF,  v  = 
ru>G ,  w 

-  H7(i/co)  ,  P  -  pB7(wi/)  ,  f  =  zj(uv)  with  the  following  boundary  conditions: 

1)  at  f  -  rj0:  F  -  0,  G  -  1,  H  -  a,  B  -  0 

2)  at  f  -  r,0-s5m;  F  -  C  -  0 


reduced  to  the  following  system  of  equations  [33] . 
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Where  r)Q  is  the  axial  position  of  the  interface  and  5ro  is  the  momentum 
boundary  layer  thickness,  a  is  the  suction  parameter  as  defined  earlier.  The 
solution  to  the  above  equations  depends  on  the  suction  parameter,  i.e.,  for  a 
given  suction  parameter  the  flow  field  can  be  calculated  once  and  used  for 
various  species  calculations. 

To  obtain  better  numerical  resolution  in  both  radial  and  axial  direction,  the 
species  concentration,  equation  (5)  in  steady  state,  is  nondimens inal ized  by 
two  reference  lengths.  The  axial  distance  is  nondimensional ized  by  the 
crystal  radius  while  the  axial  distance  is  nondimens ional ized  by  Up/D.  The 
resulting  equation  in  steady  cylindrical  polar  coordinates  is. 
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with  boundary  condition: 

ack 

1)  — —  -<  1  Kokl-0 


at  C  -  H 


°  D 


2) 


ck  =  c,  k  (0 


at  C'  =  (  n0  -  Sc)  — 


3) 


ac^ 

ar- 


=  0 


a  ck 

4)  - —  =o 

ar- 


at  f  -  0 


at  t"  =  1 


where  w'  -  w/u^  -  H(f)/a,  u'  -  u/Rsco,  f'  -  zu'/D  “  aScf ,  r'  -  r/Rs .  (r') 

is  the  species  concentration  at  the  lower  boundary  of  the  thin  layer  region 
for  specie  k,  and  must  be  matched  with  the  values  obtained  from  "bulk"  flow 
calculations.  If  is  a  constant,  the  solution  for  reduces  to  one 
dimensional  solutions. 
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V 


3.  NUMERICAL  ANALYSIS 


3-1.  Introduction 


The  numerical  procedure  used  to  solved  the  governing  equations  is  a  consis¬ 
tently  split  linearized  block  implicit  scheme  originally  developed  by  Briley 
and  McDonald  [35]  and  embodied  in  a  computer  code  termed  MINT,  an  acronym  for 
Multidimensional  Implicit  Nonlinear  Time  -  Dependent .  The  basic  algorithm  has 
been  further  developed  and  applied  to  both  laminar  and  turbulent  flows.  Since 
the  scheme  has  been  described  in  detail  in  several  publications  available  in 
the  open  literature,  it  will  not  be  detailed  here.  Rather  only  a  brief  outline 
of  the  procedure  will  be  given  in  the  following: 

The  governing  equations  are  replaced  by  an  implicit  time  difference 
approximation,  optionally  a  backward  difference  or  Crank-Nicolson  scheme. 

Terms  involving  nonlinearities  at  the  implicit  time  level  are  linearized  by 
Taylor  series  expansion  about  the  solution  at  the  known  time  level,  and  spatial 
difference  approximations  are  introduced.  The  result  is  a  system  of 
multidimensional  coupled  (but  linear)  difference  equations  for  the  dependent 
variables  at  the  unknown  or  implicit  time  level.  To  solve  these  difference 
equations,  the  Douglas-Gunn  procedure  for  generating  alternating-direction 
implicit  (ADI)  splitting  schemes  as  perturbations  of  fundamental  implicit 
difference  schemes  is  introduced  in  its  natural  extension  to  systems  of  partial 
differential  equations.  This  ADI  splitting  technique  leads  to  systems  of 
coupled  linear  difference  equations  having  narrow  block-banded  matrix 
structures  which  can  be  solved  efficiently  by  standard  block-elimination 
methods . 

The  method  centers  around  the  use  of  a  formal  linearization  technique  adapted 
for  the  integration  of  initial -value  problems.  The  linearization  technique, 
which  requires  an  implicit  solution  procedure,  permits  the  solution  of  coupled 
nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree  of 
accuracy)  by  a  one-step  noniterative  scheme.  Since  no  iteration  is  required  to 
compute  the  solution  for  a  single  time  step,  and  since  only  moderate  effort  is 
required  for  solution  of  the  implicit  difference  equations,  the  method  is 
computationally  efficient;  this  efficiency  is  retained  for  multidimensional 
problems  by  using  ADI  matrix  splitting  techniques.  The  method  is  also 
economical  in  terms  of  computer  storage,  in  its  present  form  requiring  only  two 
time  levels  of  storage  for  each  dependent  variable.  Furthermore,  the  splitting 
technique  reduces  multidimensional  problems  to  sequences  of  calculations  which 
are  one -dimensional  in  the  sense  that  easily- solved  narrow  block-banded 
matrices  associated  with  one-dimensional  rows  of  grid  points  are  produced. 
Consequently,  only  these  one -dimensional  problems  required  rapid  access  storage 
at  any  given  stage  of  the  solution  procedure,  and  the  remaining  flow  variables 
can  be  saved  on  auxiliary  storage  devices  if  desired.  Since  each  one 
dimensional  split  of  the  matrix  produces  a  consistent  approximation  to  the 
original  system  of  partial  differential  equations,  the  scheme  is  termed  a 
consistently  split  linearized  block  implicit  scheme.  Consistent  splitting  has 
been  shown  by  a  number  of  authors  [36]  to  considerably  simplify  the  application 
of  the  intermediate  split  boundary  conditions. 

The  present  work  involved  an  extension  of  the  LBI  scheme  :  in  every  time  step, 
a  boundary  fitted  coordinate  system  is  constructed  to  accommodate  the  time 
dependent  characteristic  of  the  free  surface  and  melt-crystal  interface.  A 


simple  Eulerian  height  function  is  used  to  represent  the  free  surface  and  the 
melt-crystal  interface.  This  function  is  defined  as  the  distance  from  a 
reference  line  and  is  obtained  from  the  kinematic  conditions  equation  (12)  and 
equation  (14).  The  computational  domain  is  shown  in  figure  3.1.1. 

For  numerical  calculations,  the  height  function  is  discretized  into  points 
sometime  called  markers.  These  markers  are  then  used  to  define  the  time 
dependent  boundaries.  Since  a  boundary  fitted  coordinate  system  is  used,  the 
boundary  is  part  of  the  grid  system  at  all  times.  For  this  reason,  flow 
variables  at  the  moving  boundary  are  well  defined  and  application  of  the 
boundary  conditions  at  these  boundaries  is  straightforward. 

\  The  numerical  procedure  can  be  summarized  in  four  steps:  1)  the  governing 

•  equations  of  each  phase  and  their  corresponding  boundary  conditions  are 

|  transformed  into  a  boundary  fitted  coordinate  system  according  to  the  shapes 

of  the  interfaces;  2)  the  equations  of  each  phase  are  solved  in  the 
transformed  coordinates  using  the  LBI  scheme;  3)  interfacial  positions 
(include  free  interface  and  melt-crystal  interface)  are  advanced  using  the 
i  kinematic  conditions,  and  4)  advance  time  step  and  repeat  the  sequence.  A 

flow  chart  for  the  present  numerical  scheme  is  shown  in  figure  3.1.2. 

3-2.  Transient  Capability 

The  solution  algorithm  for  iterative  solution  of  steady  solutions  has  been 
i  written  as  the  following  split  LBI  scheme  (equations  A- 8) 
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which  is  an  ADI  scheme  to  approximate  the  linearized  unsplit  scheme  (equation 
A-6) . 

(A  -  PAtLn)  <«>n  +  1  -  «>")  =  At  (Dn  +  s n)  (35) 


to  second  order  in  At.  Combining  equations  (34a)  through  (34c)  gives 
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A  comparison  of  equations  (35)  and  (36)  shows  that  splitting  error  terms  such 
as 

(pAt)2  L1a“1l2  (4>n+1  -  <t>n)  (37) 


arise  due  to  the  ADI  scheme.  Although  these  splitting  error  terms  vanish  for 


a  steady  solution  (since  ^n+1  -  <j>n)  ,  they  do  affect  transient  accuracy  in 
unsteady  applications.  The  term  (equation  (37))  is  of  second  order  in  At 
and,  hence,  does  not  lower  the  formal  accuracy  from  that  of  equation  (35). 
However,  the  error  term  (equation  (37))  has  a  different  functional  form  from 
the  truncation  error  of  the  unsplit  scheme,  and  when  matrices  such  as  A'1 
are  stiff  (i.e.,  widely  differing  eigenvalues),  the  splitting  error  can  cause 
large  errors  in  transient  solutions.  Stiffness  in  the  matrix  A"1  occurs 
at  near  incompressible  Mach  numbers. 

The  splitting  error  can  be  removed  by  iterating  at  each  time  step.  One  method 
of  iterating  is  to  replace  ^n+1  by  +  &4> ^  in  the  unsplit  scheme  (35)  ,  divide 
by  At,  and  add  an  artificial  time  derivative  approximation  A  A 4>^/^T 
to  equation  (35)  Here,  k  is  the  iteration  index,  A r  is  an  artificial  time 
step,  and  A  may  be  chosen  to  be  different  from  A.  The  result  can  be  written 
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A  Douglas-Gunn  type  of  splitting' is  introduced,  and  the  result  can  be  written 
as 


(B  -  Pl")  A<t>  *  ■=  RHS  of  Eq.  (38) 
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where  B  ~  An/Ar  +  An/At. 

Note  that  when  the  iteration  equations  (39a)  through  (39c)  converges  to  a 
steady  solution  (A^^O)  then  ,  and  the  steady  solution  satisfies 

the  linearized  unsplit  scheme  (equation  (35)).  The  steady  solution  is 
independent  of  A  and  A r,  and  these  quantities  can  be  chosen  to  obtain  rapid 
convergence.  The  matrix  conditioning  technique  of  [36]  can  be  employed  for 
low  Mach  number  calculations,  and  then  A  is  replaced  by  A  multiplied  by  a 
diagonal  conditioning  matrix. 

Since  an  iteration  is  being  performed  at  each  time  step,  there  is  the  option 
of  updating  the  time  linerizations  which  were  used  in  equation  (35).  In  this 
case,  equation  (38)  is  replaced  by 
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and  split  in  the  same  manner  as  equations  (39a)  through  (39c). 


If  equation  (40)  is  used,  then  the  converged  solution  will  satisfy  the 
nonlinear  implicit  scheme 
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Test  calculations  have  been  performed  using  equations  (39a)  through  (39c)  on 
simple  problems,  and  it  was  found  that  convergence  at  each  time  step  can  be 
obtained  in  approximately  3  to  15  iterations,  depending  on  the  physical  time 
step  used.  Fewer  iterations  are  required  as  the  physical  time  step  At  is 
reduced . 

3 - 3 .  Body-Fitted  Coordinate  Transformation 

The  numerical  solution  of  the  Navier-Stokes  equations  for  problems  with  a 
moving  boundary  is  complicated  by  the  need  to  simultaneously  compute  both  the 
flow  field  and  the  time-dependent  boundary  location  as  the  calculation 
proceeds  in  time.  In  calculations  of  this  type,  stress  conditions  or  heat 
flux  conditions  are  applied  along  the  time  dependent  boundaries.  A 
time -dependent  body-fitted  coordinate  system  is  therefore  most  suitable  and 
chosen  for  the  present  studies,  because  stress  conditions  and  heat  flux 
conditions  can  be  specified  along  the  coordinate  lines. 

Body- fitted  coordinate  systems  avoid  the  need  for  irregular  shaped  finite 
difference  computational  cells  at  computational  boundaries  thereby  alleviating 
the  complex  computer  coding  problems  associated  with  irregular  shaped  cells. 

In  addition,  spatial  truncation  error  associated  with  irregular  finite 
difference  cells  is  avoided  by  the  use  of  body-fitted  coordinates. 

Various  approaches  are  available  for  generating  body- fitted  coordinates, 
although  in  general  they  can  be  categorized  into  two  generic  groups  - 
algebraic  and  elliptic  method.  Algebraic  methods  have  the  advantage  of  being 
conceptually  simple  for  complex  boundary  geometries.  However,  it  is  often 
necessary  to  establish  suitable  functional  forms  for  each  family  of 
coordinates  for  each  new  configuration.  The  main  advantage  of  elliptic 
generation  is  that  the  method  is  applicable  to  complex  geometries  without  the 
need  to  modify  the  basic  approach  for  each  new  configuration. 

In  the  present  work,  an  algebraic  method  is  used  for  all  calculations  except 
the  case  with  a  mensicus.  In  the  meniscus  case,  elliptic  methods  based  on  the 
work  of  Thompson  and  his  coworkers  [39,40]  are  used.  We  discuss  both  methods 
below : 

In  algebraic  method,  the  time -dependent  non-orthogonal  coordinate 
transformation  is  given  by 


where  f,  x  and  £  are  computational  coordinates  in  the  r  and  z 
directions,  respectively.  f(0 .  h(x)  are  grid  packing  transformation 
developed  by  Levy  and  Gibeling  [41]  which  have  the  form  of  hyperbolic  tangent 
function.  And 


g(C.t)-  P0  +4"  1  [Pi  ertc(-^“  (^'V)  '  (  1  (oj ))] 

^  i-1  “i 


is  a  grid  distribution  generator  originated  by  Oh  [42],  in  which  the  comple¬ 
mentary  error  function  is  used  and  0  ,  j3 .  ,  7.,  a.  are  grid  parameters  whose 
values  can  be  made  time  dependent. 

In  the  present  study,  the  grid  parameters  are  constants  in  time  and  are  chosen 
such  that  grids  are  packed  at  the  regions  where  the  flow  (and/or  temperature) 
field  changes  most  rapidly.  Since  F(0 ,  H(*)  are  constant  in  time,  the  grid 
spacings  in  the  r  and  6  direction  do  not  change  in  time  despite  the  motion  of 
the  free  surface  and  interface.  The  grid  spacings  in  the  axial  direction 
depend  on  the  instantaneous  locations  of  the  free  surface  and  melt-crystal 
interface . 


The  above  way  to  constrain  grids  is  simple,  especially  in  the  present  appli¬ 
cation  where  the  slope  of  the  free  surface  and  interface  never  becomes 
infinite.  However,  in  the  case  of  a  menicus,  the  contact  angle  is  usually 
close  to  90°  (i.e.,  the  free  surface  near  the  crystal  is  close  to  being 
vertical).  In  this  case,  the  algebraic  method  mentioned  above  will  produce 
collapsing  r  and  z  coordinate  lines  which  is  not  very  desirable.  To  avoid 
this,  we  use  elliptic  methods  [43]  to  generate  grids  for  the  case  with 
mens icus . 


3 -4 .  CRAY  Vectorization 


It  is  important  to  realize  that  any  application  program  will  probably  execute 
faster  and  be  able  to  solve  larger  problems  when  transferred  to  a  super 
computer,  e.g.,  CRAY-1S,  CRAY  X-MP.  This  is  due  to  the  fact  that  super 
computers  technology  makes  use  of  short  cycle  times  (in  the  nanosecond  range) 
and  large  memories  (in  the  Mword  range).  However,  in  many  applications,  an 
order  of  magnitude  reduction  in  execution  time  may  be  obtained  bv  paying 
attention  to  the  vectorization  of  the  code. 

The  extra  speed-up  is  obtained  by  optimally  utilizing  the  vector  processing 
capabilities  of  the  machine.  In  today's  class  VI  super  computers,  areas  which 
have  a  direct  impact  on  the  performance  of  algorithms  are  as  follows: 

a)  Vector  functional  units  =  the  vector  functional  units  in  super  computer 
improve  upon  this  performance  by  using  pipeline  vector  arithmetic.  In 
order  to  utilize  the  arithmetic  pipelines,  data  has  to  be  organized  so 
that  a  steady  stream  of  operands  can  be  supplied.  Organizing  the 
operands  into  vectors  is  one  way  to  do  this.  The  longer  the  vectors  are. 
the  more  efficient  the  pipelines  can  be  utilized. 

b)  Interface  of  vector  units  with  memory  -  the  pipelined  functional  units 
which  perform  the  vector  arithmetic,  require  a  continuous  stream  of 
operands  to  arrive  at  the  pipeline  in  order  to  achieve  optimal 
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performance.  From  Che  programmer  point  of  view,  it  is  advantageous  to 
reorganize  (if  possible)  a  computational  method  involving  two  dimensional 
arrays  so  that  the  array  is  processed  in  a  column-wise  fashion. 


The  numerical  code  MINT  has  been  under  vector ization  during  this  Phase  II 
program  and  partial  vectorization  was  achieved.  In  the  partly  vectorized 
code,  (the  version  that  was  used  through  out  these  studies),  the  run  time  was 
improved  by  12.8%  over  the  unvectorized  version. 


4.  CALCULATIONS 


4-1.  Introduction 


Considerable  effort  during  the  Phase  II  study  was  devoted  to  code  development 
for  purposes  of  extending  the  state-of-the-art  CZ  studies  to  include 
self -consistently  determined  shape  of  the  melt-crystal  interface.  In  the 
interest  of  efficient  use  of  computer  resources,  initial  test  runs  for  these 
calculations  were  performed  with  values  of  dimensionless  parameters  chosen 
primarily  for  computational  resources.  Calculations  with  parameters  relevant 
to  silicon  are  discussed  in  Sections  5  and  6. 

a)  Preliminary  Test  of  Free  Surface  Algorithms 

Preliminary  calculations  were  performed  to  check  of  free  surface  algorithm  for 
the  rotating  fluid.  The  configuration  is  the  same  as  shown  in  Figure  3.1.1 
except  the  crystal  is  absent.  In  this  calculation,  dimensional  parameters  are 
normalized  by  the  crucible  radius  and  it's  rotating  velocity.  Analytical 
solutions  with  linearized  free  surface  boundary  conditions  are  available  [44] 
for  qualitative  comparison. 

Initially,  the  free  surface  is  planar  and  the  fluid  is  quiescent  with  the 
crucible  rotating  along  the  axis  of  symmetry.  The  number  of  grid  points  used 
was  26  x  11  and  are  packed  near  the  wall  and  free  surface  (as  shown  in  Figure 
4.1.1a). 

Grids  become  nonorthogonal  when  the  free  surface  deforms  because  coordinate 
lines  are  constructed  every  time  step  according  to  the  current  shape  of  the 
free  surface. 

A  steady  state  solution  was  obtained  after  80  time  steps.  As  shown  in  Figure 
4.1.1  the  shape  of  the  free  surface  is  parabolic.  Linear  theory  [44]  provides, 

V  =  Qr  ;  u  =  w  =  0 
11  =T1o  +  \  Fr  f2 

where  u  and  w  are  radial  and  axial  velocities  respectively.  The  tangential 
zero  shear  stress  condition  is  satisfied  only  if  the  shape  of  the  free  surface 
is  nearly  planar,  i.e.  small  Fr.  In  our  calculations,  which  are  fully 
nonlinear,  Fr  -  0.5,  and  w  -  1.  Consequently,  discrepancies  on  the  shape  of 
the  free  surface  and  azimuthel  velocity  near  the  free  surface  are 
anticipated.  Maximum  differences  of  3%  for  the  free  surface  location  was 
found  at  r  -  0.6.  The  contour  lines  for  v  (Figure  4.1.1b)  near  the  free 
surface  are  perpendicular  indicating  zero  tangential  stress. 

A  second  preliminary  calculation  was  performed  to  check  the  algorithm  for  the 
case  of  axisymmetric  flow  in  a  rotating  crucible  with  a  counter-rotating 
crystal.  This  calculation  is  considered  important  because  the  algorithm  was 
substantially  generalized  after  the  Phase  I  effort,  see  e.g.,  Liu  et.al  [37]. 

The  calculations  were  performed  on  a  grid  containing  26  radial  points  and  14 
axial  points.  The  solution  converged  to  a  steady  solution  within  80  time 
steps.  While  the  Reynolds  number  used  is  100,  and  is  considered  low  in  Cz 
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processes,  it  serves  the  purpose  of  verifying  the  code.  Figure  4.1.2a  shows 
the  configuration  of  the  coordinate.  Grids  points  are  packed  horizontally  at 
regions  near  the  free  surface  and  the  rotating  crystal.  Grid  points  are 
packed  vertically  at  regions  below  the  junction  of  the  free  surface  and 
rotating  crystal,  and  near  the  crucible  wall.  Numerical  calculations  of  the 
velocity  distribution  (or  field)  are  shown  in  Figure  4.1.2b  and  Figure  4.1.2c 
and  show  good  qualitative  agreement  with  Liu,  et.al.,  [37]. 

The  computer  code  used  in  the  present  study  is  based  on  the  linearized  block- 
implicit  solution  procedure  (LBI)  for  compressible  Navier-Stokes  equations 
developed  by  Briley  and  McDonald  [35].  In  the  present  development,  the 
temperature  in  the  solid  phase,  which  is  coupled  to  the  fluid  through  the 
interfacial  shape  and  heat  flux,  is  solved  for  simultaneously  with  the  fluid 
conservation  equations.  For  this  reason,  a  stand-alone  LBI  scheme  for  the 
transient  heat  conduction  equation  in  cylindical  planar  coordinate  with  moving 
boundary  at  the  axial  direction  was  developed. 

b)  Preliminary  Test  of  Transient  Conduction  Solver 

The  next  two  examples  are  devoted  to  test  the  stand-alone  LBI  solution 
procedures.  The  first  example  is  that  of  a  semi -infinite  cylinder  immersed  in 
a  hot  bath.  The  boundary  conditions  are  shown  in  4.1.3a.  For  this  problem, 
analytical  solution  exists  and  are  obtained  from  Carslaw  and  Jaeger  [45]  as 

“  2  Jo  ( 8n  r  ) 

0S  =  1  -  2  X  exp  {  -  pn  t  ) - 11 n._ 

n=1  PnMPn) 

where  /3n,  n  -  1,2  are  the  roots  of  Bessel's  function. 


Jo(Pn)  =  0 


(Note,  do  not  confuse  y9n  with  the  volumetric  expansion  coefficient.) 

Figure  4.1.3b  shows  the  comparison  between  the  numerical  and  analytical 
results . 

Excellent  agreement  is  obtained.  The  number  of  grid  points  used  was  21  x  10 
and  At  is  0.2.  The  CPU  on  CRAY-1  is  0.22  msec  per  time  step  per  mesh 
point.  In  this  calculation,  steady  state,  as  well  as  transient  solutions  for 
the  conduction  equation  are  obtained  and  were  verified  with  analytical 
solution . 


The  next  calculation  performed  was  to  verify  the  stand-alone  LBI  procedure  wi: 
moving  boundary  conditions.  The  boundary  and  initial  conditions  are  shown  in 
Figure  4.1.4a.  Initially,  the  temperature  field  6S  and  moving  surface  at 
X  -  -1  are  obtained  from  analytical  solution.  In  every  time  step  the  moving 
surface  which  is  obtained  from  analytical  solution  is  used  as  the  boundary 
location  for  the  stand-alone  LBI  procedure.  The  numerical  solutions  can  be 
compared  to  solutions  to  the  one  dimensional  analytical  transient  heat 
conduction  equation  and  solidification  surface,  as  given  by  Carslaw  and  Jaeger 
[«]. 
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The  computed  temperature  field  is  then  compared  with  the  analytical  solution 
every  time  step.  Figure  4.1.4b  shows  the  difference  between  the  numerical  and 
analytical  solutions  at  t  -  2.2.  The  maximum  error  is  2  x  10  4  occurring  at 
the  middle  of  computational  domain.  The  number  of  grid  points  is  5  x  42  and 
At  used  is  0.05. 


4-2 .  Group  I  -  Studies  of  Melt  Dynamics 

The  purpose  of  this  group  of  calculations  is  to  study  the  hydrodynamics  of  the 
melt  under  various  rotational  conditions  and  with  different  crystal  immersion 
depths.  The  response  of  free  surface  under  these  conditions  is  also 
investigated.  All  dimensionless  numbers  in  this  group  of  calculations  were 
chosen  for  computational  purposes,  they  are  not  relevant  to  silicon.  Further, 
in  all  Phase  II  studies  the  crystal  immersion  depth,  So,  is  prespecified, 
and  fixed  in  time  and  space.  This  approximation  simplifies  the  modeling  at 
the  junction  of  the  melt-solid  interface  and  at  the  melt-free  surface 
interface.  Isothermal  conditions  are  assumed  for  the  solid  and  fluid  phases 
are  assumed  in  this  group  of  calculations,  i.e.,  the  energy  conservation 
equations  for  both  phases  are  not  solved  Three  examples,  with  different 
rotational  conditions,  u,  and  immersed  length,  So,  are  chosen  as 
representative  in  this  group.  In  this  group  of  calculations,  the  Reynolds 
number  and  Froude  number  are  chosen  ‘~o  be  100  and  0.5,  respectively.  Table 
4.2.1  shows  the  difference  of  parameters  chosen  for  the  three  examples. 

i)  Run  Tl-1 

The  first  example  is  for  a  rotating  crystal  immersed  in  a  stationary  crucible, 
i.e.  u  -  o.  Recall,  the  dimensionless  crystal  rotation  is  -1.  The  crystal 
is  immersed  in  the  melt  with  So  -  0.75.  Grid  points  are  packed  into  the 
region  below  the  free  surface  and  near  the  crystal  corner.  The  number  of  grid 
points  are  26  x  23.  Figure  4.2.1  shows  results  at  steady  state;  a) 
coordinates  configuration,  b)  secondary  flow  pattern,  and  c)  azimuthal 
velocity. 

The  shape  of  the  free  surface  is  convex  up  which  qualitatively  agrees  with  the 
linear  solution  for  rotating  cylinder.  Forced  convection  is  created  due  to 
rotation  of  crystal  and  is  shown  in  Figure  4.2.1b.  The  azimuthal  velocity 
shown  in  Figure  4.2.1c  indicates  diffusion  of  momentum  which  is  generated  by 
the  rotating  crystal. 

ii)  Run  Tl-2 

This  run  is  performed  for  the  parameters  of  Tl-1  with  the  addition  of  a 
counter  rotating  crucible.  Here  u>  -  -0.68  which  means  the  crucible  is 
counter  rotating  0.68  times  slower  than  the  crystal.  With  this  parameter 
change,  the  forced  convection  phenomena  and  shape  of  the  free  surface  are  also 
changed.  The  shape  of  the  free  surface  is  convex  up  near  the  crystal  due  to 
rotation  of  crystal  and  is  concave  up  near  the  crucible  wall  due  to  rotation 
of  crucible.  Again,  the  shape  of  free  surface  is  in  good  qualitative 
agreement  with  the  linear  solution  for  counter  rotating  concentric  cylinders. 
Vortices  are  generated  due  to  rotating  surfaces.  Note  that  the  positions  of 
the  centers  of  the  vortices  depend  on  the  relative  (is  it  relative  or  the 
actual  magnitudes)  magnitudes  of  the  rotations,  i.e.,  w.  Azimuthal  velocity 


is  shown  in  Figure  4.2.2c.  Also,  note  that  there  is  a  region  within  the  me  it 
separation  two  counter  rotating  fluids. 

iii)  Run  Tl-3 

This  run  is  performed  for  the  parameters  of  Tl-2  with  one  change,  the  crystal 
is  immersed  further  into  the  melt,  $o  —  0.25.  A  similar  shape  of  the  free 
surface  is  found.  Location  of  vortices  are  changed  due  to  geometric 
differences.  Note,  with  the  crystal  further  immersed  into  the  melt,  more 
surface  area  is  available  to  generate  shear  momentum  at  the  crystal  radius. 
This  can  be  seen  from  the  plot  of  azimuthal  flow.  Contour  lines  are  pushed 
further  towards  the  crucible  wall.  In  other  words,  the  effects  of  crystal 
rotation  is  larger  in  this  case  than  in  the  previous  case. 

For  all  the  calculations  with  free  surface,  normal  velocities  along  the  free 
surface  are  monitored  in  every  time  step.  They  are  represented  by  their 
averages  value  and  R.M.S.  value  over  the  free  surface  in  every  time  step  and 
are  plotted  in  Figure  4.2.4.  The  vanishing  of  the  R.M.S.  value  of  the  normal 
velocities  serves  as  a  criteria  of  reaching  steady  state  solution. 

With  regard  to  computational  efficiency,  for  this  group  of  calculations 
(solving  the  continuity  and  Navier-Stokes  equations  of  the  fluid  phase)  the 
CPU  on  a  CRAY-1  is  2.6  msec  per  time  step  per  mesh  point. 

4-3.  Group  II  -  Studies  of  Melt  Thermohvdraulic 

In  this  group  of  calculations,  energy  conservation  together  with  the  Navier- 
Stokes  and  continuity  for  the  fluid  phase  are  solved.  Again,  dimensionless 
numbers  are  chosen  with  respect  to  computational  conditions.  A  total  of  eight 
examples  are  performed.  The  energy  equation  for  the  crystal  phase  is  not 
solved  in  this  group  of  studies. 

Five  examples  are  performed  with  planar  free  surface  and  planar  melt-crystal 
interface.  They  were  performed  to  develop  an  understanding  of  the  phenomena 
caused  by  temperature  driven  flow.  Effects  due  to  buoyancy,  Prandlt  number 
and  heat  radiation  boundary  conditions  are  studied  in  these  five 
calculations.  The  shape  of  the  melt-crystal  interface  and  its  associated 
phenomena  are  studied  in  the  next  two  calculations,  and  the  interaction  with 
the  solid  phase  forms  the  last  study  of  this  group.  Table  4.3.1.  shows  the 
parameters  used  in  this  group  of  calculations. 

i)  Runs  T2-1 ,  T2-2,  T2-3 

Examples  with  different  Grashof  numbers  (or  nondimens ional  thermal  expansion 
coefficients)  are  used  to  study  the  transition  of  forced  convection  to  free 
convection.  And  for  reference  comparison  is  made  to  run  Tl-2.  In  these 
calculations,  crystal  and  crucible  are  counter-rotating  with  w  -  -0.68  and 
the  height  of  the  crystal  is  60  -  -0.75.  With  the  Boussineq's 
approximation  for  the  density,  the  set  of  conservation  equations  are  coupled 
through  the  gravity  force  term.  The  Prandtl  number  used  for  these  runs  is 
0.054.  Run  T2-1  to  T2-3  indicate  cases  from  forced  convection  to  free 
convection,  as  determined  by  the  choice  of  f)  equation  (7),  which  are 
respectively  0,  10  and  25  (which  corresponds  to  Gr  -  0,  1.0  x  104  ,  2.6  x 
104).  Run  Tl-2  is  identical  to  T2-1,  except  that  energy  equation  is  also 


solved  In  the  latter  case .  Identical  results  for  consequences  of  the  fact 
that  the  energy  equation  is  decoupled  from  the  momentum  equation  (fi0) . 

(Note:  In  solving  this  problem,  the  energy  momentum  equations  are  treated  as 
coupled  equations).  In  Run  T2-2,  p  is  finite  and  energy  and  moments 
equation  are  coupled,  mixed  free- forced  convection  is  found.  The  fluid  has  a 
stronger  vertical  component  near  the  wall.  In  Run  T2-3,  the  free  convection 
is  stronger  than  the  forced  convection.  The  vortex  which  was  found  under  the 
crystal  in  last  runs  disappears.  The  fluid  under  the  crystal  reverses 
direction.  When  the  forced  convection  dominates  fluid  is  pumped  out  by  the 
centrifical  force;  when  free  convection  dominates,  fluid  is  moved  toward  the 
center  by  buoyancy  forces. 

For  cases  with  rotating  crystal  in  a  stationary  crucible,  free  convection  is 
dominant  if  [14] 


In  our  case,  with  counter  rotating  crystal  and  crucible,  a  similar  criteria  is 
found  to  be 
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Note  that  the  temperature  fields  are  also  affected  by  changes  of  heat  transfer 
mechanism.  The  differences  are  not  large  because  of  the  choice  of  Prandtl 
number.  For  small  Prandtl  number,  the  heat  transfer  mechanism  is  dominated  by 
heat  conduction. 

ii)  Run  T2-4 

The  purpose  of  this  run  is  to  study  the  effects  of  Prandlt  number.  This  run 
is  identical  to  T2-2,  except  that  Prandlt  number  (see  equation  (4))  is 
increased  to  1.0.  In  run  T2-2,  heat  transfer  is  mainly  due  to  conduction  as 
can  be  seen  from  the  uniformity  of  the  temperature  contours.  Figure  4.3.4c 
shows  the  temperature  contours  of  Run  T2-4  are  lower  than  those  in  Run  T2-2, 
because  convection  effects  play  a  more  important  role  for  case  with  larger 
Prandlt  number. 

iii)  Run  T2-5 

Heat  radiation  is  imposed  at  the  free  surface  (H-10)  to  study  the  temperature 
effects  on  the  fluid.  A  constant  ambient  temperature  of  1.026  is  assumed 
above  the  free  surface.  Runs  T2-5  and  T2-4  differ  in  that  the  latter  permits 
heat  radiation  at  the  free  surface.  It  is  found  that  the  temperature  below 
the  crystal  is  lower  than  that  in  run  T2-4  (see  figure  4.3.5).  The  difference 
is  direct:  the  position  of  temperature  contours  under  the  free  surface  vary 
according  to  the  choice  of  ambient  temperature. 

iv)  Run  T2-6 

Run  T2-6  is  performed  to  study  the  effects  due  to  free  surface  conditions  with 
the  same  conditions  as  in  Run  T2-5.  These  conditions  are:  a)  counter  rotating 
crystal  and  crucible,  w  -  -0.68,  b)  mixed  forced-free  convections,  6  -  10, 


c)  heat  radiation  boundary  condition  at  the  free  surface,  H-10,  d)  large 
Prandlt  number,  Pr  -  1.0.  For  comparison  of  the  free  surface  the  results  of 
run  Tl-2  are  also  useful. 


The  shape  of  the  free  surface  found  in  this  run  is  concave  down  (see  figure 
4.3.6)  near  the  crystal  surface,  while  it  is  convex  up  in  run  Tl-2.  This  is 
because  fluid  with  lower  temperature  is  carried  downward  near  the  crystal 
surface.  The  general  shapes  of  the  free  surface  near  the  crucible  are  the 
same  in  these  two  runs.  By  allowing  the  free  surface  to  deform,  the  maximum 
magnitude  of  secondary  flow  is  larger  than  the  case  in  T2-5.  Consequently, 
more  heat  is  being  carried  down  below  the  crystal.  The  result  is  that  the 
temperature  below  the  melt-crystal  interface  is  lower  than  those  in  run  T2-5. 
Differences  in  temperature  in  this  region  will  influence  the  local  formation 
of  melt-crystal  interface,  as  will  be  shown  in  the  next  calculation.  Below 
the  crystal,  azimuthal  velocities  from  run  T2-5  and  present  run  are  also 
different,  due  to  differences  of  magnitudes  of  secondary  flow.  Figure  4.3.7 
shows  the  R.M.S.  values  of  the  normal  velocities  of  the  free  surface  as  a 
function  of  time.  Steady  state  is  reached  when  the  normal  velocities  of  the 
free  surface  are  zero. 

v)  Run  T2 -  7 

In  this  calculation  the  temperature  field  of  the  crystal  phase  is  assumed 
uniform  and  at  melt  temperature  as  in  T2-1  to  T2-6.  However,  here  the 
melt-crystal  interface  is  allowed  to  deform  according  to  the  kinematic 
condition.  This  assumption  is  valid  if  the  temperature  gradient  normal  the 
interface  is  small.  All  other  conditions  are  identical  to  the  last  run 
(T2-6).  The  dimensionless  parameters  for  this  calculation  are  S  =  22,  Pe  - 
0.1,  K  -  1  .  Indeed,  depending  upon  the  choice  of  parameters,  the  shape  of 
the  melt-solid  interface  will  deform  into  concave  or  convex.  Figure  4.3.8 
shows  the  results  at  steady  state.  The  differences  in  flow  field  and 
temperature  field  indicate  the  effect  due  to  the  deformation  of  the  melt-solid 
interface . 

This  calculation  is  the  first  to  show  the  distortion  of  the  melt  crystal 
interface  due  to  the  thermal-hydralics  of  the  melt.  There  are  two  distinct 
vorticity  patterns  as  revealed  by  the  secondary  flow  calculations.  The  fluid 
is  being  drawn  to  the  center  of  the  crystal  and  then  thrown  out  radially.  The 
fluid  is  transporting  the  heat  away  from  the  axis  of  symmetry.  There  is  a 
counter  transport  of  heat  arising  from  the  counter-rotating  crucible.  The  dip 
in  the  crystal  interface  occurs  near  the  highest  temperature  zone  of  the 
melt.  Also  note  from  the  azimuthal  distribution  that  more  of  the  melt  is 
rotating  with  the  crystal,  than  that  experienced  with  the  imposed  flat 
interface . 

4-4 .  Group  III  -  Studies  of  Species  Transport  Near  the  Interface 

As  discussed  in  Section  2-8  and  repeated  here  that  the  behavior  of  the  fluid 
adjacent  to  the  growing  crystal  interface  is  a  primary  factor  in  controlling 
the  species  uptake  to  the  crystal.  In  our  model,  we  assume  that  species 
diffusion  and  forced  convection  (due  to  rotation  of  the  crystal)  are  important 
in  a  thin  region  with  thickness  A  near  the  melt -crystal  interface.  Solute 
impurities  in  the  melt  and  from  tfie  crucible  wall  are  transported  by  "bulk" 
flow  and  reach  the  boundary  of  this  thin  region.  The  governing  equations  (2) 
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and  (3),  can  be  simplified  to  a  set  of  ordinary  differential  equations  (33]  to 
provide  a  flow  field  for  the  species  equation  (5)  in  this  region. 

The  solution  to  the  set  of  coupled  ordinary  differential  equations  equation 
(29)  through  (32)  is  a  function  of  suction  parameter,  a,  and  the  boundary 
conditions  at  f  ■  f2,  where  f  is  the  nondimens ional  thickness  of  the  layer. 
Theoretically,  in  the  case  of  a  finite  diameter  disk  rotating  in  an  infinite 
fluid,  the  radial  and  azimuthal  velocities  will  decay  to  zero  as  f  ■  ®. 
However,  in  practice,  as  shown  in  Schlichting  [3]  (with  a-o) ,  the  radial  and 
azimuthal  velocities  decay  asymptotically  to  zero  on  distance  of  within 
47(i//w) ,  i.e.  ,  |c0  -  f2|  -4. 

There  are  two  approaches  to  handle  the  bottom  velocity  boundary  conditions:  In 
the  first  approach,  we  extend  the  thin  layer  region  to  a  fictitious  length 
(i.e.,  f  -  f2  -  some  large  value),  and  apply  zero  velocity  condition 
there.  This  condition  was  shown  to  be  consistent  with  the  velocities  of  the 
bulk  flow  calculation  obtained  from  the  set  of  ordinary  differential  equations 
at  the  point  f t ,  where  f2<fi<f0-  1°  the  second  approach,  we  find  Cj  F(f,). 

G^j)  from  the  bulk  flow  calculations  and  use  them  as  boundary  conditions 
for  the  set  of  ordinary  differential  equations. 

The  first  approach  is  used  in  the  Case  1  study.  The  fictitious  length  is 
arbitrarily  taken  to  be  5  (>4) .  The  zero  condition  i.e.,  F(-5)  -  G(-5)  -  0  is 
applied  at  f2.  The  solution  to  the  set  of  ordinary  differential  equations 
in  this  case  is  shown  in  Figure  4.4.1.  The  second  approach  is  used  in  the 
Case  2  study.  The  thickness  of  the  thin  layer  |f0  -  fj|  is  obtained 
from  bulk  flow  calculations,  and  is  obtained  when  the  radial  velocity  exhibits 
a  maximum  as  it  is  searched  axially  downward  from  the  interface.  F(f1), 

G(fj)  are  then  obtained  by  averaging  u  and  v  in  the  radial  direction, 
respectively.  It  is  found  that  -  0.83  and  F(fj)  -  0.16,  G(fj)  -  0.58. 
Results  with  this  condition  are  shown  in  Figure  4.4.2.  The  suction  parameter 
is  taken  to  be  zero  in  both  cases.  By  comparing  the  results  from  Cases  1  and 
2,  we  find  that  the  solution  to  the  set  of  ordinary  differential  equation 
using  the  first  approach  matches  the  solution  from  the  bulk  flow  calculations 
for  fj<f<f0,  where  fj«-l  and  valid  at  the  approximation  used. 

Physically,  it  means  that  the  fluid  characteristic  within  J(u/u>)  away  from 
the  melt-crystal  interface  can  be  described  using  the  assumption  of  the  thin 
layer  model . 

By  assuming  that  the  amount  of  solute  impurities  is  small,  i.e.,  the  solute 
redistribution  dot'  not  affect  melt  hydrodynamics,  the  convective-diffusive 
species  equation  can  be  decoupled  from  the  melt  hydrodynamics  equation.  In 
the  following  analysis  for  the  species  distribution,  we  will  first  use  one 
dimensional  approximation  to  compare  our  results  with  those  obtained  by  pre¬ 
vious  investigators.  To  predict  species  radial  impurities,  a  two  dimensional 
equation  together  with  bulk  flow  solution  are  used. 

a)  One -Dimens ional  Species  Distribution 

The  one -dimens ional  governing  equation  for  the  species  distribution  is 
obtained  by  simplying  equation  (33),  i.e., 


with  boundary  conditions,  9C 

'(1  k°,C  =  0  31  ;'=;o 

C=CL  at  C'  =  C'c 

The  nondimens ional  velocity  w'  is  obtained  from  the  solution  of  two  dimensional 
axisymmetric  flow  equations  (29)  through  (32).  Here,  f  '0  is  the  location  of 

the  interface  and  can  conveniently  be  taken  to  be  0,  and  f 'c  is  equal  to  aScfc, 
where  fj<f^<5"0.  Within  the  range  of  parameters  (a,  k0 ,  Sc)  for  silicon 
crystal  growth,  we  found  that  C  approaches  Cl  if  ir^l>  a^c •  In  the  following 
calculations,  |f^,|  is  taken  to  be  2aSc-  a  longer  computational  domain 
than  necessary.  This  is  chosen  to  ensure  that  the  downstream  boundary  condi¬ 
tion  do  not  influence  the  upstream  condition  (in  this  case,  C  at  the 
interface) .  There  are  various  materials  (with  different  kG)  in  the  melt 
and  calculations  were  performed  to  cover  the  range  of  kQ  from  0.0004(In)  to 
1.4  (0).  Figure  4.4.3  shows  the  solute  distribution  in  the  melt  for  various 
k0  using  the  axial  flow  field  obtained  from  Case  1.  Figure  4.4.4  shows  the 
solute  distribution  in  the  melt  for  various  k0  using  the  axial  flow  field 
obtained  from  Case  2.  In  this  case,  -  0.83  aSc. 

From  the  crystal  grower's  point-of-view,  solute  impurities  distribution  in  the 
melt  is  of  less  interest.  It  is  more  desirable  that  the  solute  impurities 
level  in  the  crystal  are  known.  Here  a  commonly  used  parameter  is  the 
effective  distribution  coefficients,  keff,  k0  C0/Cl,  where  kQ  is 
the  equilibrium  distribution  coefficient  defined  as: 


p  E  concentration  of  the  species  in  the  crystal  at  the  interface 
concentration  of  the  species  in  the  melt  at  the  interface 

k0  is  usually  taken  from  the  phase  diagram.  The  keff  is  then  a  measure  of 
concentration  level  of  a  particular  specie  in  the  crystal  if  the  concentration 
level  in  the  melt  is  known.  According  to  Burton,  et  al.  (BPS)  [28],  keff 
can  be  written  as, 

k  =  _ _ 

e  <  +  (  1  -  k0  )  exp  ( -  A  ) 

where  a  =  -  In  [  1  aSc|”exp  (  Sc  H  cK)  d$] 

o  o 

=  -  In  [  1  -  aSc  F(a,  Sc  )  ] 

BPS'  expression  for  kef£  is  exact  with  respect  to  the  one  dimensional  solute 
conservation  equation  and  boundary  conditions.  The  function  F(a,Sc)  depends 
on  the  hydrodynamics  of  the  melt  (i.e.,  H(f'))  and  can  be  obtained  by 


numerical  methods.  Table  4.4.1  shows  numerical  correlations  for  4  obtained 
by  previous  investigators.  Their  relations  are  limited  to  small  values  of 
aSc^/3  or  (aScV^).  For  a  particular  value  of  a  and  Sc,  9  different 
values  of  k0  (ranging  from  0.0004  to  1.4)  are  used  to  determine  the  solute 
concentrations  at  the  interface,  i.e.,  keff  are  known  (keff  -  k0 
C0/Cl) .  From  BPS'  expression  for  keff,  a  set  of  A  (for  each  kc)  is 
obtained.  As  are  then  averaged  and  a  data  point  is  generated  for  that 
particular  a  and  Sc.  Within  the  range  of  parameters  for  silicon  crystal 
growth,  we  used  various  suction  parameters,  a,  (from  0.0005  to  0.1)  and 
Schmidt  numbers,  Sc,  (from  10  to  54)  to  determine  C0  and  consequently 
keff.  From  the  numerical  data  we  are  able  to  find  a  numerical  expression 
for  F(a , Sc ) ; 

A  =  -  In  {  1  -  aSc  ( -j-  +  2Sc  0  375 exp  (-  1.52  aSc0  55)  )  ] 

This  formula  can  be  reduced  to  the  BPS  correlation  for  small  aSc^/^.  Figure 
4.4.5  shows  the  comparison  of  our  empirical  relations  with  others.  The 
comparison  is  good. 

b .  Two  Dimensional  Species  Distribution 

In  the  last  section,  we  assumed  that  the  radial  variation  of  species  distri¬ 
bution  is  small  compared  to  the  axial  variation.  In  this  section  we  do  not 
make  this  assumption  i.e.,  the  two  dimensional  axisymmetric  species  equation 
(33)  is  used  for  the  numerical  calculations.  The  radial  and  axial  velocities 
of  the  melt  are  obtained  from  the  solution  of  the  reduced  Navier- Stokes 
equations  (29)  through  (32).  The  solute  distribution  level  (r')  at  f ' 

-  aSc,  which  is  the  boundary  condition,  is  needed  and  obtained  from  the 
solution  of  the  'bulk'  flow  calculations  equation  (5).  Two  calculations  are 
performed  to  validate  the  computer  code. 

The  first  calculations  assumed  Cl(t')  at  f'  -  aSc  be  a  constant.  The 
purpose  of  this  calculations  is  to  show  that  calculations  using  a  two  dimen¬ 
sional  equation  should  reproduce  the  results  of  the  calculations  using  the  one 
dimensional  equation  if  the  axial  boundary  conditions  do  not  vary  in  the 
radial  direction.  The  result  is  shown  in  Figure  4.4.6.  In  this  case  we  used 
a  -  0.05,  Sc  -  53  and  kQ  -  0.3  x  5(P).  Physically  it  simulates  the 
distribution  of  dopant  phosphorous  which  has  an  equilibrium  distribution 
coefficient  of  0.35  in  the  melt  near  the  melt-crystal  interface.  The  crystal 
is  being  grown  in  a  rate  of  lmm/min  and  rotating  at  4rpm.  At  the  melt-crystal 
interface,  the  normalized  concentration  level  is  about  1.78.  The  specie 
distribution  (in  the  axial  direction)  is  identical  to  the  result  obtained  from 
one  dimensional  calculations. 


In  the  second  case,  we  allow  at  f'  -  aSc  varies  in  the  radial 
direction.  The  values  of  Cl  at  f '  -  aSc  are  obtained  from  the  'bulk' 
flow  calculations  equation  (5).  The  radial  and  axial  velocities  of  the  melt 
are  obtained  from  the  solution  of  the  reduced  Navier- Stokes  equations  (29) 
through  (32).  The  purpose  of  this  calculation  is  to  demonstrate  the 
capability  of  the  thin  layer  model  in  predicting  the  solute  distribution  (in 
both  axial  and  radial  directions)  in  the  melt  near  the  melt-crystal 


interface.  Nondimens ional  parameters  used  are  the  same  as  in  first  case.  The 
result  is  shown  in  Figure  4.4.7.  It  is  found  that  the  nonuniform  species 
distribution  caused  by  natural  convection  has  an  effect  on  the  distribution  of 
the  species  at  the  melt-crystal  interface. 
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5.  CALCULATIONS  II 


5 - 1 .  Group  IV  -  Applications  to  Silicon  Crystal  Growth 

In  the  previous  section  we  verified  the  numerical  procedures  which  have 
implemented  in  the  existing  computer  code.  The  calculations  were  performed 
using  nondimens ional  parameters  not  applicable  to  silicon  crystal  growth. 

These  parameters  were  chosen  to  minimize  the  computation  resources  needed  to 
verify  the  code.  In  this  section,  we  summarize  the  numerical  calculations 
performed  with  parameters  relevant  to  Cz  silicon  growth  under  normal 
conditions.  Since  not  all  growth  parameters  are  measurable,  e.g.  heat 
transfer  coefficient  at  the  free  surface,  etc.,  approximations  were  used  for 
those  parameters. 

A  total  of  six  two  dimensional  (no  9  dependence)  calculations  were  performed 
for  silicon  CZ  growth.  The  first  four  cases  are  simulations  of  silicon 
crystal  growth  with  a  crystal  diameter  of  8.8cm  in  a  13cm  diameter  crucible. 
The  silicon  is  rotating  at  1.32rpm,  while  the  crucible  is  counter-rotating  at 
0.5rpm.  The  height  of  the  grown  crystal  is  taken  as  5.1cm  and  the  height  of 
the  melt  (between  the  crystal  interface  and  crucible  bottom)  is  taken  to  be 
3.8cm.  The  first  two  cases  assumed  the  crucible  has  a  uniform  temperature  of 
1462  C,  while  the  third  and  fourth  cases  assumed  a  crucible  temperature  of 
1432  C.  The  fifth  case  shown  is  a  simulation  of  silicon  crystal  growth  in  a 
26.5cm  diameter  crucible.  The  crystal  diameter  is  the  same  as  in  the  previous 
cases,  and  the  crucible  is  counter-rotating  at  a  rate  of  0  13rpm.  This  was 
chosen  to  give  the  same  crucible  Reynolds  number  as  for  the  previous  case. 

The  height  of  the  melt  is  taken  to  be  17.7cm  in  this  case,  and  the  crystal 
interface  is  assumed  to  be  at  the  equilibrium  Si  melting  temperature  1412  C. 
The  heat  radiation  condition  at  both  the  melt-argon  and  crystal-argon 
interfaces  is  applied  in  all  the  cases.  The  sixth  case  simulated  a  growth  of 
a  large  diameter  crystal  (28.58cm).  In  this  case  the  Reynold's  number  is 
about  7  times  larger  than  the  other  cases. 

The  purpose  of  Cases  1  and  2  is  to  compare  the  effects  of  the  temperature 
variation  along  the  melt-crystal  interface  with  and  without  the  assumption  of 
a  planar  interface.  Case  3  demonstrates  that  a  near  planar  interface  (or 
nearly  uniform  temperature  variation  along  the  interface)  can  be  achieved  by 
reducing  the  bouyancy  effects,  and  proper  adjustment  of  the  heat  radiation 
condition.  Case  4  shows  the  capability  to  calculate  a  case  with  a  meniscus. 
Further  details  of  these  calculations  are  presented  below.  Case  5  is  to  study 
the  effects  of  large  crucible,  while  Case  6  shows  a  growth  of  large  diameter 
crystal . 

In  all  the  cases  studied,  a  common  flow  pattern  is  observed:  at  least  two 
vortexes  exist  in  the  melt.  One  vortex  is  found  near  the  crystal.  This 
vortex  has  a  radially  outward  flow  due  to  the  rotation  of  the  crystal.  The 
radial  velocity  (near  the  crystal  interface)  increases  as  the  radius  of  the 
crystal  and  reaches  a  maximum  at  about  one  crystal  radius.  It  diminishes  as 
it  goes  further  outward.  At  some  radial  distance  (usually  greater  than  one 
crystal  radius)  the  flow  turns  downward,  reversed  its  direction  and  flow 
radial  inward.  This  change  of  direction  is  caused  by  another  vortex  which  is 
driven  by  the  natural  convection.  The  point  where  the  first  vortex  turns 
direction  depends  on  the  relative  strength  of  the  forced  convection  (induced 
by  the  rotation  of  the  crystal)  and  the  natural  convection  (induced  by  the  hot 


crucible  walls) .  The  strength  of  the  forced  convection  depends  on  rotational 
rate  and  radius  of  the  crystal.  While,  the  strength  of  the  natural  convection 
depends  on  the  magnitude  of  the  Grashof  number  and  Prandlt  number.  The 
Grashof  number  depends  on  the  physical  properties  of  the  melt,  temperature 
difference  between  the  crucible  wall  and  the  melt  temperature,  and  a 
characteristic  length.  The  Prandlt  number  depends  on  the  physical  properties 
of  the  melt.  For  this  reason,  the  flow  pattern  in  a  large  crucible  (Case  5) 
is  much  more  complicated. 

Details  of  each  case  study  follows: 

CASE  1 

In  this  calculation,  the  melt-crystal  interface  is  assumed  planar,  i.e.,  the 
shape  of  the  melt-crystal  interface  is  predetermined.  The  assumption  is  valid 
if  the  net  heat  flux  along  the  interface  is  approximately  constant.  The 
average  temperature  of  the  argon  gas  (above  the  free  surface)  is  assumed  to  be 
1355  C.  Emissivities  for  the  crystal  surface  and  free  surface  are  taken  to  be 
0.2  (a  constant  for  simplicity).  Figure  5.1.1  shows  the  grid  used  for  the 
numerical  simulations.  Grid  points  are  packed  near  the  interface  and  crystal 
corner  to  obtain  better  resolution  in  those  regions. 

Since  the  melt-crystal  interface  is  planar,  the  rotating  crystal  acts  as  a 
centrifugal  disk  pump  which  sucks  up  the  melt  and  spins  it  out  in  the  radial 
direction.  As  can  be  seen  from  Figure  5.1.2,  an  Ekman  layer  of  thickness 
O.llcra  is  formed  under  the  planar  interface.  The  theoretical  value  from  a 
finite  diameter  rotating  disk  in  a  semi -infinite  fluid  without  free  convection 
is  about  0.14cra.  The  radial  outward  flow,  due  to  the  rotating  crystal,  meets 
the  radial  inward  flow,  due  to  natural  convection  driven  by  the  hot  crucible 
wall,  at  a  radial  distance  1.2cm  away  from  the  crystal  corner.  Temperature 
contours  for  both  phases  are  also  shown  in  Figure  5.1.2.  The  zero  heat  flux 
condition  has  been  used  at  the  top  of  the  crystal,  hence  a  certain  axial 
distance  (in  this  case  5.1cm),  the  temperature  within  the  crystal  becomes 
uniform.  The  axial  heat  flux  along  the  melt-crystal  interface  is  nonuniform, 
which  is  a  consequence  of  the  neglect  of  the  latent  heat  of  crystallization  at 
the  interface. 

Figure  5.1.3  shows  the  heat  fluxes  along  the  interface  in  the  melt  and  crystal 
phases.  The  difference  between  the  curves  indicates  the  net  heat  flux  removed 
from  the  heat  gained  due  to  latent  heat  of  crystallization.  Due  to  natural 
convection,  the  higher  temperature  melt  is  carried  towards  the  interface  for 
r>0.3.  This  motion  of  the  melt  produces  nonuniform  heat  flux  along  the 
interface  in  the  melt  phase  and,  consequently,  nonuniform  solidification  at 
the  interface  results.  Physically,  nonuniform  solidification  along  the 
interface  means  a  nonplanar  interfacial  shape.  This  example  shows  that  if  the 
assumption  of  a  planar  interface  is  used,  uniformity  of  the  net  heat  flux  to 
the  interface  must  be  checked  to  verify  the  consistency  of  the  assumption. 

The  concentration  distribution  of  an  arbitrary  impurity  is  shown  in  Figure 
5.1.4.  The  distribution  is  calculated  by  assuming  no  loss  or  creation  of 
species  at  any  surface.  For  convenience,  the  concentration  level  is 
normalized  to  a  reference  concentration.  Since  the  species  boundary  layer 
thickness  is  about  Sc  -  7.3  times  thinner  than  the  momentum  boundary  layer 
thickness,  numerical  resolution  of  the  concentration  distribution  will  be 


inadequate  near  the  interface.  A  zero  condition  is  used  at  the  interface  for 
this  reason,  hence,  in  this  'bulk'  flow  calculation,  species  concentration 
levels  near  the  interface  are  valid  only  outside  the  species  boundary  layer. 

A  submodel  (thin  layer)  was  developed  to  resolve  the  species  boundary  layer 
near  the  melt-crystal  interface  as  discussed  in  Section  4-4.  Concentration 
levels  at  the  interface  from  the  'bulk'  flow  calculation  are  used  as  boundary 
conditions  in  the  submodel  calculation.  In  the  thin  layer,  natural  convection 
effects  are  assumed  to  be  small,  and  the  fluid  flow  is  solely  due  to  the 
rotation  of  the  crystal.  In  this  case,  the  thickness  of  the  thin  reg'.on  is 
taken  to  be  0.8 J(i//u>)  (0.11cm).  Phosphorus,  which  has  a  segregation 
coefficient,  k0 ,  of  0.35,  is  chosen  for  the  numerical  calculations.  Using 
the  rotational  rate  of  the  crystal  and  a  calculated  pulling  rate  equation  (14) 
of  0.2mm/min,  the  suction  parameter  'a.'  is  determined  to  be  0.0175. 

Figure  5.1.5  shows  the  concentration  distribution  inside  the  thin  layer.  We 
note  again  that  from  a  crystal  grower's  point  of  view,  the  species 
concentration  distribution  in  the  melt  is  of  no  interest;  rather,  the  species 
concentration  level  in  the  crystal  is  of  primary  interest.  The  commonly  used 
effective  distribution  coefficient  keff,  defined  as  k0*C0/CL,  is  a 
measure  of  concentration  level  of  a  particular  specie  in  the  crystal.  Here, 
k0  is  the  equilibrium  distribution  coefficient  which  is  a  ratio  of  specie 
concentration  in  the  crystal  at  the  interface  to  the  concentration  of  specie 
in  the  melt  at  the  interface,  and  k0  is  usually  taken  from  the  phase 
diagram.  For  a  particular  specie,  if  the  concentration  of  the  specie  in  the 
melt  near  the  interface  is  known,  keff  can  be  determined.  From  our 
numerical  calculations,  the  keff  (based  on  C0  and  Cl)  for  phosphorus 
varies  from  0.427  to  0.408,  while  a  numerical  correlation  using  one 
dimensional  approximation  by  BPS  [34]  gave  a  value  of  0.445.  The  one 
dimensional  analysis  [34]  assumed  a  quadratic  variations  for  the  axial 
velocity  and  ignored  the  radial  variation  of  the  species  concentrations. 

CASE  2 

The  geometric  configuration  in  this  calculation  is  the  same  as  in  the  first 
case,  except  that  the  shape  of  the  interface  is  also  determined  as  part  of  the 
solution  from  the  heat  flux  condition  at  the  melt-crystal  interface.  The 
crystal  pulling  velocity  can  be  determined  from  the  area  average  net  heat 
transfer  across  the  interface,  since  this  is  a  necessary  condition  for  a 
steady  state  solution  to  exist.  In  other  words,  for  the  physical  parameters 
used,  the  crystal  is  growing  at  a  constant  diameter.  The  interfacial  position 
is  considered  as  steady  state  if  the  local  pulling  velocity  along  the 
interface  is  within  3%  of  the  average  pulling  velocity  obtained  from  average 
net  heat  available  for  solidification. 

Figure  5.1.6  shows  the  coordinates  at  the  end  of  the  steady  state 
calculation.  The  initial  coordinate  mesh  is  orthogonal,  which  is  the  same  as 
that  shown  in  Figure  5.1.1,  but  due  to  deformation  of  the  interface  the 
coordinates  become  nonorthogonal .  Figure  5.1.7  shows  the  secondary  flow  and 
temperature  contours  where  it  is  seen  that  the  vortex  structure  and 
temperature  distributions  are  similar  to  those  for  Case  1.  The  shape  of  the 
interface  becomes  convex  for  the  parameters  used,  and  the  maximum  int  rface 
deformation  is  about  0.9cm.  The  pulling  velocity  obtained  from  the 
calculations  is  about  0.2mm/min.  Local  heat  fluxes  in  the  melt  and  crystal 


phases  along  the  melt-crystal  interface  are  plotted  in  Figure  5.1.8.  The 
difference  between  the  heat  fluxes  indicates  the  amount  of  heat  gain  from 
crystallizaticn. 

In  this  example,  we  have  shown  that  nonuniform  solidification  at  the  interface 
leads  to  nonplanar  interface  and  nonuniform  temperature  near  the  melt-crystal 
interface . 


CASE  3 

The  purpose  of  this  calculation  is  to  show  that  a  'near'  planar  interface  can 
be  achieved  by  adjusting  various  parameters.  In  this  particular  calculation, 
the  crucible  temperature  has  been  reduced  such  that  the  strength  of  natural 
convection  is  reduced  to  provide  a  more  uniform  temperature  distribution  near 
the  melt-crystal  interface.  On  the  other  hand,  rotational  rate  of  the  crystal 
is  increased  to  provide  a  stronger  forced  convection  with  a  similar  result. 

In  this  calculation,  the  temperature  of  the  crucible  is  assumed  to  be  1432  C. 
Emissivities  for  the  free  surface  and  crystal  surface  are  taken  to  be  0.2  and 
0.08,  respectively.  As  seen  from  Figure  5.1.9,  the  influence  of  natural 
convection  is  substantially  decreased.  Temperature  contours  for  both  phases 
are  also  shown.  The  maximum  deviation  of  interfacial  height  is  less  than 
1.2%,  while  the  maximum  deviation  of  the  free  surface  height  (without 
considering  surface  tension)  is  less  than  0.2%.  Heat  fluxes  of  both  phases 
along  the  interface  are  shown  in  Figure  5.1.10,  where  again  the  difference 
between  the  two  curves  indicates  the  amount  of  heat  gained  from 
crystallization.  Figure  5.1.11  shows  the  normalized  net  heat  fluxes  along  the 
interface  with  and  without  the  assumption  of  planar  interface,  where  it  is 
seen  that  the  heat  flux  is  much  more  uniform  in  the  case  with  the  calculated 
nonplanar  interface. 

The  specie  concentration  distribution  of  an  arbitrary  impurity  is  shown  in 
Figure  5.1.12.  The  boundary  conditions  used  are  the  same  as  in  Case  1.  The 
more  uniform  distribution  near  the  interface  which  is  found  in  this  case 
confirms  that  the  effect  of  natural  convection  in  transporting  species  from 
the  melt  toward  the  interface  is  reduced  with  the  parameters  chosen  for  this 
simulation. 

Again,  a  thin  layer  submodel  calculation  near  the  melt-crystal  interface  is 
performed  to  resolve  the  species  boundary  layer  (see  Figure  5.1.13).  Using 
the  rotational  rate  of  the  crystal,  and  the  calculated  crystal  pulling 
velocity  (0 . 04mm/min) ,  the  suction  parameter  is  determined  to  be  0.0035.  From 
the  present  numerical  calculations,  the  kegf  for  phosphorous  only  varies 
from  0.369  to  0.368,  due  to  the  nearly  uniform  concentration  distribution 
while  the  one  dimensional  approximation  by  BPS  [34]  is  0.368. 

CASE  4 

The  purpose  of  this  calculation  is  to  demonstrate  the  numerical  capability  of 
the  existing  code  to  handle  a  case  with  a  meniscus.  In  this  calculation  a 
Bond  number  of  58  and  contact  angle  of  15  are  used.  The  Bond  number  and  the 
contact  angle  are  determined  from  the  physical  properties  of  silicon.  The 


free  surface  shape  is  determined  by  the  surface  tensional  force  and  inertial 
force  due  to  rotations.  Other  parameters  are  the  same  as  in  Case  3  except  the 
following:  the  average  temperature  of  the  argon  gas  is  assumed  to  be  1379  C. 

Due  to  the  sharp  change  in  slope  of  the  free  surface  near  the  crystal,  a  new 
method  of  constructing  the  grid  is  used,  with  the  coordinates  shown  in  figure 
5.1.14.  The  boundary- fitted  coordinates  are  generated  using  elliptic  partial 
differential  equations  based  on  the  work  of  Thompson,  et  al.,  [39,40].  The 
origins  of  the  elliptic  generation  technique  are  based  on  an  analogy  between 
the  requirements  of  two  dimensional  coordinate  systems,  and  the  properties  of 
streamline  and  potential  function  distributions.  This  grid  generation 
technique  have  been  modified  and  applied  successfully  [43]  in  complicated 
geometries.  The  temperature  contours  of  both  phases  and  secondary  flow  are 
shown  in  Figure  5.1.15.  The  maximum  deviation  of  the  interfacial  height,  in 
this  case  is  within  4%  (i.e.,  ~lmm)  while  the  maximum  deviation  of  the  free 
surface  is  about  7mm.  The  calculated  meniscus  shape  when  the  rotational 
effects  are  not  considered  have  been  compared  with  the  analytical  solution  of 
Young- Laplace ' s  equation  [20].  The  comparison  is  in  good  agreement.  The  flow 
pattern  obtained  in  this  case  is  quite  similar  from  those  in  previous  cases. 
The  vortice  which  under  the  free  surface  in  the  previous  case  is  now  'pushed' 
into  the  interior  of  the  melt  under  the  crystal  by  the  meniscus. 

The  general  flow  pattern,  however,  does  not  change:  an  Ekman  layer  under  the 
crystal;  radially  outward  flow  caused  by  the  rotating  crystal  is  forced  to 
turn-around  by  the  radially  inward  flow  caused  by  natural  convection;  axial 
flow  beneath  the  crystal  center  is  upward  or  downward  depends  on  whether  it  is 
above  or  below  the  Ekman  layer. 

CASE  5 

The  calculation  shown  in  (Cases  1  through  4)  are  done  in  a  crucible  of 
relatively  small  diameter.  Note  that  the  short  height  of  the  crucible  means 
that  the  simulations  are  done  for  the  near-end  stage  of  a  nonreplenishable 
process.  In  this  case,  a  larger  size  crucible  is  used.  The  diameter  of  the 
crucible  is  26.5cm  with  a  melt  height  of  17.7cm  and  the  length  of  the  crystal 
is  assumed  to  be  13.3cm.  The  rotational  rate  of  the  crystal  is  the  same  as 
before,  but  the  rotational  rate  of  the  crucible  is  reduced  so  that  the 
crucible  Reynold's  number  is  the  same  as  for  the  previous  cases.  In  physical 
terms,  the  crucible  is  rotating  at  0.13rpm,  and  the  temperature  of  the 
crucible  is  assumed  to  be  1432  C. 

Figure  5.1.16  shows  the  coordinates  for  this  calculation  where  there  are  66 
radial  and  60  axial  mesh  points  in  the  melt  phase,  and  42  radial  and  22  axial 
points  in  the  crystal  phase.  Figure  5.1.17  shows  the  secondary  flow  and 
temperature  distributions,  while  a  detailed  velocity  vector  plot  near  the 
melt-crystal  interface  is  shown  in  Figure  5.1.18.  In  this  calculation, 
several  more  vortexes  are  found  in  the  melt  fluid.  The  most  distinguished  one 
is  the  one  centered  at  6.8cm  from  the  crucible  bottom.  It  carries  the 
hot-melt  fluid  to  a  region  near  the  crystal  interface.  Several  recirculating 
vortexes  are  also  observed.  As  shown  in  Figure  5.1.18,  an  Ekman  layer  of 
approximately  0.11cm  is  also  found.  We  observed  the  similar  results  as  in  the 
previous  cases  since  the  nondimensional  number  in  the  governing  equations  are 
the  same.  However,  there  is  no  geometric  similarities  with  the  last  cases 


(e.g.,  aspect  ratio  of  the  crucible  are  not  the  same  as  before).  In  this 
calculation,  the  ambient  temperature  of  the  argon  gas  is  adjusted  so  that  a 
melt-crystal  interface  of  a  maximum  derivation  of  about  1%  is  obtained. 

CASE  6 

The  calculation  shown  in  this  case  is  to  simulate  the  crystal  growth  process 
to  yield  a  large  diameter  crystal.  A  crystal  of  diameter  19.05cm  is  grown 
crucible  of  diameter  28.58cm.  The  height  of  the  melt  is  9.53cm.  The  crystal 
is  rotating  at  2rpm  while  the  crucible  is  counter-rotating  at  0.5rpm.  In 
nondimensional  terms,  the  Reynolds  number,  Res,  based  on  the  crystal  radius 
and  rotation  is  7000  and  the  Reynold's  number  Rec,  based  on  the  crucible 
radius  and  rotation  is  4000. 

Figure  5.1.19  shows  the  coordinates  for  this  calculation  where  there  are  68 
radial  and  48  axial  mesh  points  in  the  melt  phase,  and  35  radial  and  25  axial 
points  in  the  crystal  phase.  Figure  5.1.20  shows  the  secondary  flow  and 
temperature  contours  of  the  melt  and  crystal  phases.  In  this  calculation,  the 
forced  convection  induced  by  both  the  crystal  and  crucible  rotations  is  larger 
than  the  previous  cases  and  this  results  in  a  much  more  complicated  flow 
pattern  and  temperature  distribution  in  the  melt.  The  radial  velocities  near 
the  crucible  are  much  larger  than  in  the  previous  cases.  These  large  radial 
velocities  push  the  counter-rotating  vortex  upward  and  closer  to  the  crucible 
side  wall.  The  center  of  the  vortex  is  at  about  one  crystal  radius  in 
contrast  to  about  0.9  radius  in  previous  cases.  The  strong  radial  flow  makes 
a  strong  return  flow  at  the  free  surface  near  the  crucible  side  wall.  This 
strong  return  flow  brings  the  hot  melt  from  the  bottom  and  side  crucible  wall 
to  the  free  surface . 
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6.  THREE  DIMENSIONAL  CALCULATIONS 


6- 1 .  Group  V  -  Applications  to  Continuous  Feed  Silicon  Growth 

A  simulation  of  Czochralski  silicon  growth  with  a  Continuous  Liquid  Feed  (CLF) 
option  has  used  to  demonstrate  the  three  dimensional  capability  of  the 
computer  code  developed.  A  continuous  CLF  crystal  growth  was  designed  to 
permit  crystal  growth  from  a  small  melt  volume.  It  has  been  reported  that 
large  melt  volume  changes  during  growth  of  the  crystal  was  a  strong  influence 
to  the  crystal  quality.  The  advantage  of  small  melt  volume  used  in  CLF  is 
then  obvious.  Another  advantage  of  the  CLF  method  is  the  high  crystal  yield 
per  crucible  i.e.,  longer  crystal  can  be  grown  from  each  crucible. 

A  description  of  the  set-up  of  the  CLF  is  following.  Two  identical 
Czochralski  crystal  growth  chamber  are  used.  One  crucible  (the  melter)  serves 
as  a  supplier  of  melt  fluid  to  the  other  (the  grower).  These  two  chambers  are 
isolated  from  each  other  to  eliminate  problems  associated  with  SiO  build-up  in 
the  growth  chamber.  All  growth  conditions  in  these  two  chambers  are  identical 
except  that  the  vertical  motion  of  the  two  solid  rods  are  in  opposite 
directions.  The  melt  fluid  in  the  two  crucibles  is  connected  by  an  inverted 
U-tube  at  the  free  surface  and  melt  fluid  is  transferred  from  the  melter  to 
the  grower  by  hydrostatic  pressure  differentials  through  the  inverted  U-tube. 

Since  we  are  interested  in  the  physics  inside  the  grower  chamber,  numerical 
calculations  are  done  for  the  growth  chamber  only  and  the  three  dimensional 
aspect  of  the  problem  is  obtained  by  assuming  a  specific  temperature 
distributions  within  a  small  region.  The  justification  for  this  temperature 
assumption  is  as  follows.  The  contact  surface  between  the  connecting  tube  and 
free  surface  is  modeled  as  an  area  with  a  temperature  different  from  that  of 
the  free  surface.  The  flow  rate  in  the  inverted  U-tube  is  assumed  to  be  small 
compared  to  the  melt-fluid  velocity  (secondary  velocity)  near  the  free 
surface.  The  temperature  of  the  contact  surface  is  specified  as  following. 

In  CLF  growth,  the  environment  in  the  grower  chamber  is  the  same  as  the 
environment  in  the  melter  chamber,  except  that  the  directions  of  the  growth 
are  reversed  in  the  two  chambers.  The  temperature  profile  at  the  inlet  of  the 
inverted  U-tube  (in  the  melter  side)  should  approximately  be  equal  to  the 
temperature  profile  at  an  area  on  the  free  surface  of  a  grower  crucible.  Due 
to  the  passing  of  melt  the  fluid  inside  the  tube,  the  outlet  temperature  of 
the  tube  can  approximate  as  the  average  of  the  inlet  temperature  profile. 

In  our  numerical  simulation,  an  rectangular  area  of  0.87  x  0.84cm2  is 
assumed  to  be  the  contact  surface  of  the  connecting  tube.  This  contact 
surface  is  located  about  1cm  from  the  crucible  wall  at  the  free  surface.  The 
initial  conditions  correspondence  to  a  two  dimensional  Czochralski  growth  with 
boundary  conditions  as  stated  in  Case  3,  Section  4.3.  Once  the  three 
dimensional  steady  state  solution  is  obtained,  we  change  the  temperature  of 
the  contact  area  to  the  average  of  the  temperature  in  the  area.  Calculations 
with  a  rectangular  spot  of  constant  temperature  at  the  free  surface  is  then 
performed.  Figure  6.1.1  shows  the  coordinates  in  the  r -8  plane  for  the  melt 
and  solid  phases.  Notice  that  the  grids  are  packed  where  the  connecting  tube 
is  in  contact  with  the  melt  phase  at  the  free  surface.  Coordinates  in  the  r-z 
plane  are  the  same  as  in  the  two  dimensional  cases.  Figure  6.1.2  shows  the 
velocity  plots  in  the  r -9  planes  at  various  heights.  The  free  surface  is 
located  at  about  f  —  0.84.  The  small  temperature  perturbation  due  to  the 


connecting  tube  does  not  seem  to  have  a  large  effect  on  the  azimuthal 
velocities.  Figures  6.1.3  and  6.1.4  shows  the  temperature  contours  of  the 
melt  and  crystal  phases  respectively.  The  cooler  temperature  melt  is  well 
mixed  at  0.6cm  below  the  free  surface.  Further,  the  small  temperature 
disturbance  does  not  have  influence  as  the  temperature  contours  of  the  crystal 
Figure  6.1.5  and  6.1.6  show  the  secondary  flow  and  temperature  contours  of  the 
melt  phase  under  the  contacting  area.  In  order  to  look  at  how  the  temperature 
disturbance  distribute  in  space,  we  plot  the  non-dimensional  temperature 
difference  in  the  radial  direction  for  various  height.  The  plot  is  shown  in 
figure  6.1.7.  As  shown  in  Figure  6.1.8,  the  temperature  disturbance  behaves 
like  an  exponential  decay  function  in  height. 


-  41  - 


7.  CONCLUSIONS  AND  RECOMMENDATIONS 


A  numerical  algorithm  to  follow  free  and  melt-crystal  interfaces  in 
conjunction  with  the  LBI  scheme  has  been  developed  and  tested  for  Czochralski 
crystal  growth  problems.  The  physical  domain  is  transformed  using  a  time 
dependent  boundary  fitted  coordinate  system  and  the  nonlinear  interfacial 
boundary  conditions  are  applied  along  these  boundaries  which  can  be  evolving 
in  time.  The  solid  and  liquid  phases  are  coupled  through  the  heat  flux 
conditions  along  the  melt- solid  interface.  The  velocity,  temperature,  and 
species  distribution  of  the  melt  phase,  temperature  distribution  of  the 
crystal  phase,  along  with  the  positions  of  the  melt-crystal  interface  and  free 
surface,  are  obtained  by  the  consistency  split  LBI  scheme.  The  crystal 
pulling  velocity  is  also  determined  from  the  necessary  condition  of  the 
interface  tracking  algorithm.  A  thin  layer  submodel  which  can  provide  species 
concentration  levels  near  the  melt-crystal  interface  is  also  included.  The 
procedure  has  been  utilized  for  sample  calculations  which  demonstrate  the 
feasibility  of  the  approach  and  its  potential  value  for  assessing  the 
Czochralski  crystal  growth  process. 

With  the  numerical  capability  developed  in  the  present  work,  the  thermal 
variation  and  dopant  concentration  homogeneity  along  the  melt-crystal 
interface  can  be  predicted  based  on  the  external  growth  parameters. 

Moreover, a  crystal  pulling  velocity  can  be  provided  for  uniform  growth. 

Since  our  numerical  scheme  is  written  in  primitive  variable  forms,  extension 
to  study  three  dimensional  effects  is  not  limited  and  was  demonstrated  in  a 
simulation  of  Continuous  Liquid  Feed  (CLF)  Czochralski  silicon  growth. 

In  the  present  studies,  we  found  that  with  the  physical  properties  of  the 
molten  silicon  and  crucible  sizes  used  in  the  industry,  the  shape  of  the  free 
surface  caused  by  rotation  of  the  crucible  and  crystal  cannot  be  approximated 
as  planar.  Surface  tensional  force  at  the  melt-crystal-gas  junction  and  near 
the  free  surface  is  important.  The  former  causes  the  existence  of  a  mensicus 
and  latter  causes  Marangoni  flow.  The  shape  of  the  melt-crystal  interface  is 
important  to  know  because  it  indicates  the  uniformity  of  the  heat  flux 
(directly  the  uniformity  of  the  thermal  stress)  in  the  grown  crystal.  To 
achieve  a  near  planar  melt-crystal  interface,  the  rotational  rate  of  the 
crystal  and/or  the  heat  transfer  conditions  of  the  argon  gas  need  to  be 
determined. 
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APPENDIX  A 


Split  LBI  Algorithm 

Linearization  and  Time  Differencing 

The  system  of  governing  equations  can  be  written  as  a  single  grid 
point  in  the  following  form: 

3  H(<j>)/3t  =  DU)  +  S(4>)  (A- 1 ) 

where  $  is  the  column-vector  of  dependent  variables,  H  and  S  are 
column-vector  algebraic  functions  of  and  D  is  a  column  vector  whose 
elements  are  the  spatial  differential  operators  which  generate  all  spatial 
derivatives  appearing  in  the  governing  equation  associated  with  that 
element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximations  of  (A-l): 

(Hn+1  -  Hn)/At  =  8(Dn+1  +  Sn+1)  +  (1-0)  (Dn  +  S")  (A-2) 

where,  for  example,  Hn+*  denotes  H($n+*)  and  At  *  tn+l  -tn.  The 
parameter  0  (0.5  8  1)  permits  a  variable  time-centering  of  the  scheme, 

with  a  truncation  error  of  order  (A-  ,  (0  -  1/2)  AtJ. 

A  local  time  linearization  (Taylor  expansion  about  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear 
differential  operator  L  such  that 

0n+l  -  D"  +  L"  (♦"*'  -  ♦")  ♦  0  (it2)  (A- 3) 

Similarly, 

Hn+1  =  H°  (3H/3$)"  ( $n+ '  -$")■*•  0  (At2)  (A-4) 

S"+I  -  Sn  ♦  (3S/3*)n  (*"+‘  -  ♦")  ♦  0  (At2)  (A-5) 

-  46  - 


Equations  (A-3,  4,5)  are  inserted  into  Eq.  (A-2)  to  obtain  the  following 
system  which  is  linear  in 

(A  -  8At  L°)  ($n+1  -  $")  -  At  (d°  +  S°)  (A-6) 

and  which  is  termed  a  linearized  block  implicit  (LBI)  scheme.  Here,  A 
denotes  a  square  matrix  defined  by 


A  =  (3H/3<J>)n  -  BAt  (3S/3^)n  (A-7) 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (A-6)  is 
split  using  ADI  techniques.  To  obtain  the  split  scheme,  the 
multidimensional  operator  L  is  rewritten  as  the  sum  of  three 
“one-dimensional"  sub-operators  (i  =  1,  2,  3)  each  of  which  contains 
all  terms  having  derivatives  with  respect  to  the  i-th  spatial  coordinate. 
The  split  form  of  Eq.  (A-6)  can  be  derived  either  by  following  the 
procedure  described  by  Douglas  and  Gunn  in  their  generalization  and 
unification  of  scalar  ADI  schemes,  or  using  approximate  factorization.  In 
either  case,  for  the  present  system  of  equations  the  split  algorithm  is 
given  by 

(A  -  BAclJ)  (<fr*  -  $n)  ■  At  Dn  +  Sn)  (A-8a) 

(A  -  BAtLj )  (♦**-  «|>n)  ■  A  ($*  -  <f.n)  (A-8b) 

(*  -  Baa")  {,ntl  -♦”)-*  (♦**-  ♦“)  (a-bc) 

where  and  $**  are  consistent  intermediate  solutions^.  If  spatial 
derivatives  appearing  in  and  D  are  replaced  by  three-point  difference 
formulas,  then  each  step  in  Eqs.  (A=8a),  b,  c)  can  be  solved  by  a 
block-tridiagonal  elimination. 
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Heat  transfer  model  at  crystal  surface 

If  the  heat  convection  boundary,  as  outlined  in  section  2-6,  is  used,  the 
Nusselt  number  must  be  known.  The  boundary  condition  is  applied  along  the 
exposed  crystal  surface.  The  most  suitable  heat  transfer  model  can  be  found  in 
Bird  [64]  -  Heat  transfer  coefficients  for  forced  convection  around  submerged 
objects.  A  correlation  is  available  under  the  following  assumptions. 

a)  uniform  surface  temperature 

b)  gas  properties  are  evaluated  at  film  temperature 

c)  gas  approaching  velocity  is  constant 

d)  gas  flow  is  cross  stream 

log  (Nu  Re'1  Pr'1/3)  =  m  log  Re  +b 

where  m  =  -0.4725 
b  -  -0.3147 

and  the  Nusselt  and  Reynolds  numbers  are  nondimensionalized  by  the  diameter  of 
the  cylinder. 


Group 

Description 

A 

Studies  on  hydrodynamics  of  the  melt  phase  with 

planar  melt-crystal  interface  and  free  surface 

B 

Studies  on  temperature  distribution  of  the  melt 

and/or  crystal  phase 

C 

Combined  problem  of  melt  and  crystal  phase 

Melt  Phase 

Crystal  Phase 

Authors 

Method 

Convection 

Heat 

Radiation 

Free 

Surface 

Melt 

Interface 

Method 

Heat 

Radiation 

Forced 

Free 

Kobayashi 
and  Arizumi 
1970 

FDM 

2-D 

Re<500; 

Pr=0.3 

Rotating  Crystal 

Gr=60 

£  =0.2 

Planar 

N/A 

N/A 

N/A 

Langlois 
and  Shir 
1977 

FDM 

2-D 

Re=2; 

Pr=2 

Counter  Rotation 

Gr=<3 

N/A 

Planar 

N/A 

N/A 

N/A 

Langlois 

1977 

FDM 

2-D 

Re=40; 
Pr=0.08 
Counter  Rotation 

Gr=9x105 

N/A 

Planar 

N/A 

N/A 

N/A 

Kobayashi 

1978 

FDM 

2-D 

Re=40; 
Pr=0.01 .  0.1 
Rotating  Crystal 

3  5 

10  <Gr<10 

Insulated 

Planar 

N/A 

N/A 

N/A 

Langlois 

1979 

FDM 

2-D 

Re=1700; 
Pr=0.03 
Rotating  Crystal 

Gr<4.7x10? 

£=0.33 

Planar 

N/A 

N/A 

N/A 

Kobayashi 

1980 

FDM 

2-D 

Re=40; 

Pr=1 

Counter  Rotation 

Gr=500 

Insulated 

Planar 

N/A 

N/A 

N/A 

Kobayashi 

1981 

FDM 

2-D 

40<Re<200; 

Pr=1 

Counter  Rotation 

Gr<105 

Insulated 

Planar 

N/A 

N/A 

N/A 

Crochet,  et  al. 
1983 

FEM 

2-D 

Red  000; 
Pr=0.01 

Counter  Rotation 

Gr<106 

Insulated 

Planar 

N/A 

N/A 

N/A 

Arizumi  and 
Kobayashi 
1972 

FDM 

2-D 

N/A 

N/A 

e  =0.2 

Planar 

Included 

FDM 

e  =0.2 

Williams  and 
Reusser 
1983 

N/A 

N/A 

N/A 

N/A 

N/A 

Planar 

FEM 

Specified 

Temperature 

Derby  and 
Brown 
1985 

FEM  for 
LEC 
2-D 

N/A 

N/A 

e=0.8* 

Y  =0.05 

Meniscus 

Included 

FEM 

£=0.8* 

Y  =0.05 

Ramachandrar 

1985 

N/A 

N/A 

N/A 

N/A 

N/A 

Included 

FEM 

£=0.64 

Chang  and 
Brown 
1984 

FEM  for 
Bridge  man 
2-D 

Pr=1 

No  Rotation 

5  6 

2x10  <Gr<10 

N/A 

N/A 

Included 

FEM 

N/A 

Liu,  et  al. 
1982 

FDM 

2-D 

Re=4000; 
Pr=0.08 
Counter  Rotation 

Gr=0,  2x10  2 

N/A 

Planar 

N/A 

N/A 

N/A 

Present 

Work 

FDM 

2-D  &  3-D 

100<Re<4000 
Pr=1,  0.054 
Counter  Rotation 

6 

0<Gr<6x10 

e  =0.2 

Included 

& 

Meniscus 

Included 

FDM 

£  =0.2 

*  q  *e(04  -  04)  +y(0-  0a) 


TABLE  1.1.2  -  Numerical  Studies  of  Cz  Crystal  Growth  by  Various  Authors 


RUN  # 
PARAMETERS 


TI-2 


TI-3 


U) 


Tl  -I 


0  -0.68 


0.75  0.75 


-0.68 


025 


RUN 


TABLE  4.3.1  -  Parameters  Used  in  Group  II  Studies 


TABLE  4.3.2  -  Features  Studied  in  Groups  I  and  II 


CASES  1  through  5 
Res  =  1000. 

Rec  =  850. 

Pr  =  0.054 
Fr  =  0.001 

Gr  Tm 

- =  5.0 

Re  2  AT 

Pe  =  27. 


TABLE  5.1.1.  Nond linens ional  Parameters  Used  in  Group  IV  studies 
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FIG.  1.1.1  -  Basic  Principle  of  Czochralski  Crystal  Growth  Method 
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FOR  THE  MELT  REGION 
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CONSTRUCT  COORDINATES 
FOR  THE  SOLID  REGION 
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ADVANCE  TIME  STEP 


NO  /  STEADY 
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(  STOP  J 


FIG.  3.1.2  -  Flow  Chart  of  the  Numerical  Method 


-  58  - 


w 


0.505 


FIG.  4.1.1 


Free  Surface  Calculations  in  a  Rotating  Fluid 
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DEFINITION 


Comparison  Between  the  Numerical  and  Analytical  Results 
of  a  Semi-Infinite  Cylinder  Immersed  in  a  Hot  Bath 


ERROR 


FIG.  4.1.4  -  Comparison  Between  the  Numerical  and  Analytical  Results 
of  the  One-Dimensional  Transient  Heat  Conduction  with  a 
Solidification  Line 
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FIG.  4.2.4 


Average  and  RMS  Values  Along  the  Free  Surface 
vs.  Time  in  a  Typical  Group  I  Study 
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FIG.  4.3.3  -  Run  T2-3:  P  ~  25,  H  -  0,  Pr  -  0.054 
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FIG.  4.3.4  -  Run  T2-4:  B  -  10,  H  -  0,  Pr  -  1.0 
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-  Run  T2-7 :  P  -  10,  H  -  10,  Pr  -  1.0,  with  Free 
Surface,  Melt  Interface  a'nd  Crystal  Phase 


FIG.  4.3.8 
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FIG.  4.4.1  -  Normalized  Velocities  in  the  Thin  Region  Model 


'*m 

v| 

ft' 


I 


I 


0.06 


0.00 


0.75 


0.50 


0.00 


■5.00  -3.75  -2.50  -1.25  0.00 


£  *  ijui/v 


FIG.  4.4.2  - 


Comparison  of  Two  Normalized  Velocities  Calculations: 

a)  Assume  Quiescent  Condition  at  £  -  5; 

b)  'Bulk’  Flow  Conditions  at  £  -  0.83 
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FIG.  4. 4. 3. a  -  Species  Distributions  for  Various  Dopants 
Using  Velocity  Field  from  Fig.  4. 4. 2. a 
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FIG.  4 . 4 . 3 . b  -  Species  Distributions  for  Various  Dopants 
Using  Velocity  Field  from  Fig.  4. 4. 2.  a 
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FIG.  4 . 4 . 4 . b  -  Species  Distributions  for  Various  Dopants 
Using  Velocity  Field  from  Fig.  4.4.2.b 
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Group  IV:  Case  1  -  Secondary  Flow  Field  and  Temperature  Contours 
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Group  IV:  Normalized  Species  Distribution  in  the  Thin  Region 
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FIG.  5.1.6  -  Group  IV:  Case  2  -  Coordinates 
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Group  IV:  Case  2  -  Secondary  Flow  Field  and  Temperature  Contours 
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FIG.  5.1.8  -  Group  IV:  Case  2  -  Temperature  Gradients  of 
the  Phases  Along  the  Melt-Crystal  Interface 
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FIG.  5.1.12  -  Group  IV:  Case  3  -  Normalized  Species 
Distribution  in  the  Melt  Phase 
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,1.16  -  Group  IV:  Case  5  -  Coordinates 
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FIG.  6.1.1.  -  Group  V:  3-D  -  Coordinates 
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Group  V : 3-D  -  Azimuthal  Velocities  in  Various  Horizontal  Planes 


FIG.  6.1.5  -  Group  V:  3-D  -  Secondary  Flow  Field  in  a  Vertical 
Plane  Cutting  Across  the  Perturbed  Temperature  Spot 
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FIG.  6.1.8  -  Group  V:  3-D  -  Temperature  Difference 
at  r  -  1.2  vs.  Depth  of  the  Melt 


